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Propagation  characteristics  of  an  electromagnetic 
(EM)  wave  inside  a  waveguide  are  greatly  modified  by 
the  introduction  of  dielectrics  into  the  guide. 
Frequency  of  cutoff,  attenuation,  and  power  flow 
distribution  are  all  properties  of  the  EM  wave  that 
are  highly  dependent  on  the  physical  structure  and 
composition  within  the  guide.  Shielding  applications 
take  advantage  of  the  large  amount  of  attenuation 
provided  by  a  waveguide  when  the  incident  EM  wave 
has  a  frequency  lower  than  the  cutoff  frequency  of 
the  waveguide.  In  practice,  small  sections  of 
waveguides  are  inserted  through  the  metal  walls  of  a 
shielded  enclosure  to  construct  an  air  passage  that 
does  not  compromise  the  shielding  effectiveness  of  the 
shelter.  Typically,  however,  other  fluids  and  materials 
besides  air  must  be  transferred  inside  the  enclosure. 
Therefore,  hoses  are  sometimes  inserted  through  these 
waveguides,  which  can  reduce  the  shielding  provided 
by  the  waveguide  above  and  below  cutoff. 

This  work  has  investigated  the  EM  propagation 
through  these  loaded  waveguides  structures,  placing 
special  emphasis  on  wave  attenuation.  The  study  was 
divided  into  two  problems:  the  modal  coupling  and 
excitation  of  a  finite  inhomogeneous  guide. 

Solutions  to  the  radially  inhomogeneous  infinite 
circular  waveguides  have  been  obtained  numerically 
from  an  exact  theoretical  development  beginning  with 
first  principles.  These  inhomogeneous  waveguides  are 
layered  with  circular,  concentric  annuli,  possibly 
having  complex  constitutive  parameters.  The  modes 
that  were  found  to  exist  in  these  waveguides  are 
hybrid,  meaning  that  they  have  both  axial  E-  and  H- 
Helds.  The  foundations  and  numerical  solutions 
developed  were  verified  experimentally. 
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CHAPTER  1 
INTRODUCTION 


The  propagation  characteristics  of  an  electromagnetic  (EM)  wave  inside  a 
waveguide  are  greatly  modified  by  the  introduction  of  dielectrics  into  the  guide.  The 
frequency  of  cutoff  (where  the  wave  begins  propagation),  the  attenuation  and  the 
distribution  of  power  flow  are  all  properties  of  the  EM  wave  that  are  highly  dependent  on 
the  physical  structure  and  composition  within  the  guide.  Shielding  applications  take 
advantage  of  the  large  amount  of  attenuation  provided  by  a  waveguide  when  the  incident 
EM  wave  has  a  frequency  lower  than  the  cutoff  frequency  of  the  waveguide.  In  practice, 
small  sections  of  waveguides  are  inserted  through  the  metal  walls  of  a  shielded  enclosure  to 
construct  an  air  passage  which  does  not  compromise  the  shielding  effectiveness  of  the 
shelter.  Typically,  however,  other  fluids  and  materials  besides  air  must  be  transferred 
inside  the  enclosure.  To  accomplish  this,  hoses  are  occasionally  inserted  through  these 
waveguides;  hence,  the  shielding  provided  by  the  waveguide  above  and  below  cutoff  can 
be  reduced.  It  is  the  aim  of  this  work  to  investigate  the  EM  propagation  through  these 
loaded  waveguide  structures  placing  special  emphasis  on  the  attenuation  of  the  wave. 

This  loaded  waveguide  problem  may  be  separated  into  two  main  subproblems,  (i) 
the  modal  coupling  and  excitation  of  a  finite  inhomogeneous  guide,  and  (ii)  the  modal 
solutions  and  wave  properties  for  an  infinite  inhomogeneous  guide.  These  two  parts  are 
related  by  the  fact  that  once  the  excited  modes  are  identified  in  the  finite  guide,  the 
attenuation  and  other  attributes  of  the  wave  can  be  calculated  using  the  infinite  guide 
solutions.  Identifying  the  excited  modes  is  a  difficult  task  theoretically  and  is  not 
considered  in  this  work.  Experimentally,  however,  a  large  number  of  measurements  were 
performed  on  such  an  arrangement.  Instead,  this  work  emphasizes  the  wave  properties  not 
directly  dependent  on  the  source  excitation,  but  concentrates  more  on  the  inhomogeneous 
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waveguide  penetration  itself  and  the  characteristic  properties  of  the  geometry  and  loading 
materials. 

As  a  first  approximation  to  the  theoretical  analysis  of  this  problem,  the  waveguide 
structure  is  modelled  as  an  infinite  waveguide  containing  layered  concentric  dielectrics 
forming  homogeneous  annuli.  In  this  vein,  there  has  been  much  work  done  for  the  two 
concentric  dielectric  problem  in  which  a  dielectric  rod  is  inserted  inside  the  guide. 
Pincherle  [1]  was  one  of  the  first  to  consider  this  twc-dielectnc  guide  and  found  solutions 
for  both  the  rectangular  and  circular  guides.  For  the  inter,  only  azimuthally  symmetric 
modes  were  considered.  The  two  papers  by  Clarricoats  [2]  provide  a  more  complete, 
authoritative  exposition  on  the  two -dielec trie  guide  and  they  also  build  the  stage  from  which 
many  papers  over  the  next  decade  ( 1960's)  on  the  topic  of  the  rod-insert  waveguide  were 
set.  Recently,  a  number  of  papers  concerning  the  two-dielectric  insen  waveguide  and 
resonators  have  appeared  by  Zaki  and  Aria  [3]  and  Zaki  and  Chen  [4],  [51.  In  these  papers, 
more  emphasis  is  placed  on  the  nenazimuthally  symmetric  modes  and  their  corresponding 
field  patterns  and  resonant  frequencies  for  red-loaded  cavities  of  similar  cross  sections. 
Until  recently,  very  little  if  anything  has  been  published  for  the  cases  in  which  more  than 
two  dielectrics  fill  the  waveguide.  Bruno  and  Fridges  [o]  consider  a  variant  of  the  three- 
dielectric  problem —  that  of  a  two-dielectric  lossless  rod  guide.  Chou  and  Lee  [7]  analyze 
the  case  of  a  perfect  electrically  conducting  (pec)  waveguide  with  multiple  coatings  having 
perhaps  complex  permeabilities.  This  is  related  to  the  problem  at  hand,  but  here,  however, 
lossy  dielectric  fillings  rather  than  lossy  magnetic  coatings  are  relevant.  In  this  work,  the 
theoretical  formulation  and  subsequent  numerical  solutions  for  the  waveguide  filled  with 
any  number  of  concentric,  perhaps  lossy,  dielectrics  will  be  presented.  Both  electric  and 
magnetic  losses  are  allowed  or  any  combination  of  p^s  and  e/s;  however,  only  electrical 
losses  will  be  considered  here. 

Solutions  to  this  loaded  waveguide  problem  hinge  on  finding  roots  of  a 
characteristic  equation  which  is  derived  from  first  principles.  The  roots  of  this  equation  are 
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the  axial  wavenumbers,  kj,  for  a  given  azimuthal  variation.  There  are  an  infinite  number  of 
these  roots  for  each  azimuthal  variation  and  developing  a  method  by  which  a  root  can  be 
assigned  to  a  mode  is  a  large  part  of  the  total  solution  process.  One  method  of  identifying 
roots  in  the  inhomogeneous  case,  and  one  which  is  used  extensively  in  this  work,  is  to 
"trace”  the  root  in  the  complex  plane  as  some  parameter  is  continuously  varied.  The 
starting  point  for  this  tracing  process  is  chosen  where  the  mode  type  and  root  are  known, 
usually  the  homogeneous  guide  case,  and  then  the  1^  value  is  continuously  computed  while 
some  physical  parameter  is  incremented,  such  as  the  radius  of  the  dielectrics,  the 
constitutive  parameters  of  the  dielectrics,  the  frequency,  or  perhaps  the  amount  of  loss. 
This  tracing  process  becomes  especially  important,  and  also  difficult,  when  there  is  much 
loss,  since  the  roots  are  then  located  somewhere  in  the  total  complex  plane.  In  this  work, 
many  examples  will  be  given  when  very  high  dielectric  iosses  are  present,  from  which  the 
attenuation  and  axial  power  flow,  among  other  things,  will  be  examined.  In  particular,  as 
discussed  at  the  beginning  of  this  chapter,  the  three-layered,  dielectrically-loaded 
waveguide  is  very  important  in  this  study.  An  in-depth  investigation  into  the  attenuation 
for  this  inhomogeneous  guide  will  be  performed  in  Chapter  4. 

During  this  kj  tracing  process,  for  reasons  to  be  discussed  throughout  this  thesis 
beginning  in  Section  3.3,  backward-wave  regions  may  develop  for  these  dielectrically- 
loaded  waveguides.  A  backward-wave  region  is  one  in  which  the  wavefronts  are  traveling 
in  the  opposite  direction  to  the  net  power  flow.  ("Net"  power  flow  is  a  necessary  statement 
since  it  will  also  be  shown  that  the  power  flow  can  have  differing  signs  in  different 
dielectric  annuli.)  Clarricoats  and  Waldron  [8]  were  the  first  to  predict  the  existence  of 
such  a  phenomenon  in  the  two-dielectric  circular  guide.  This  topic  was  extensively  studied 
after  that  with  Clarricoats  publishing  two  excellent  papers  on  the  topic  [9],  [10]. 
Experimental  verification  for  the  existence  of  these  backward  waves,  in  addition  to  the 
theoretical  development,  was  carried  out  by  Clarricoats  and  Birtles  [11]  with  a  resonant 
cavity  technique  and  by  Clarricoats  and  Slinn  [12]  using  a  slot  in  a  loaded  waveguide. 
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Both  techniques  showed  very  good  agreement  with  the  theoretically  expected  values.  An 
extensive  list  of  applications  for  the  rod-loaded  waveguide  backward-wave  structure  was 
discussed  by  Waldron  [13].  In  this  thesis,  the  existence  of  these  backward  waves  will  be 
verified  for  the  two-dielectric  case  and  also  intensively  studied  in  the  three-dielectric  guide. 

A  new  modal  designation  scheme  for  backward-wave  modes  will  be  introduced, 
beginning  in  Section  3.5,  that  is  more  consistent  with  other  properties  of  the  wave,  in 
addition  to  the  kz,  such  as  net  power  flow,  boundedness  of  the  wave  and  continuity  in  the 
kz  trace.  This  new  scheme  relies  on  the  introduction  of  a  small  amount  of  loss  into  the 
dielectrics  in  order  that  the  proper  kz  root  can  be  chosen.  Until  the  recent  advent  of 
sophisticated  Bessel  function  computer  subroutines  which  can  compute  values  for  complex 
arguments,  this  could  not  be  done  very  easily.  Along  with  this  new  modal  designation 
scheme  comes  a  slightly  revised  look  at  the  cutoff  frequency  for  a  mode.  The  issue  at  hand 
concerning  this  topic  is  whether  kj  equals  zero  for  cutoff,  the  phase  and  group  velocities 
are  infinite  and  zero,  respectively  at  cutoff,  none  of  these,  or  maybe  some  combination. 
The  results  of  this  study  suggest  that  the  frequency  for  which  the  group  velocity  is  zero  is  a 
more  fundamental  concept  for  cutoff  than  1^=0  and,  correspondingly,  the  phase  velocity 
becoming  infinite. 

In  addition  to  backward  waves  forming  during  the  kj  tracing  process,  regions 
having  complex  kj  in  the  form  of  [3— jcx  can  also  exist,  even  in  lossless  media!  Normally, 
this  would  not  be  expected  since  wave  attenuation  in  lossless  media  with  finite  power  flow 
would  seem  to  violate  conservation  of  energy  principles.  Clarricoats  and  Taylor  [14] 
predicted  the  existence  of  these  complex  modes  in  the  two-dielectric  waveguide.  It  is 
pointed  out  there,  and  by  Chomey  [15],  that  these  complex  modes  always  exist  in  pairs 
having  complex  conjugate  axial  wavenumbers  so  that  zero  net  power  flow  is  preserved 
across  any  transverse  plane  of  the  waveguide.  Therefore,  conservation  of  energy  is 
enforced  since  no  net  axial  power  flow  is  present.  This  does  not  imply,  however,  that 
locally  there  is  no  axial  power  flow.  It  will  be  seen  later  in  this  thesis  that  the  axial  power 
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flow  is  mostly  nonzero  in  the  cross  section  of  the  guide  and  only  the  net  power  flow  goes 
to  zero.  Chomey  also  states  that  not  only  is  the  real  power  flow  zero,  but  the  complex 
modes  also  carry  no  net  reactive  power.  Laxpati  and  Mittra  [16]  also  note  the  zero  axial 
power  flow  condition  for  the  complex  modes  in  lossless  media  and  further  claim  this  to  be 
a  condition  for  these  modes  to  physically  exist  in  the  waveguide  and  not  merely  as 
"spurious,"  nonphysical  modes.  Recently,  Kalmyk,  Rayevskiy  and  Ygryumov  [17]  have 
presented  experimental  evidence  showing  the  existence  of  these  complex  modes  in  a 
circular  waveguide. 

Although  no  net  power  flow  is  carried  by  the  complex  modes  in  lossless  media, 
once  losses  are  introduced  into  the  dielectrics,  the  relative  phase  difference  between  the  two 
coupled,  complex  modes  will  make  net  power  flow  possible.  From  a  shielding  application 
standpoint,  the  amount  of  attenuation  in  'hese  complex  modes  becomes  an  important  issue. 

The  Final  topic  considered  in  this  thesis  is  the  experimental  measurements  for  the 
inhomogeneous  circular  waveguide.  There  were  two  types  of  measurement  arrangements 
used  here  which  applied  either  to  loading  dielectrics  which  were  nearly  lossless  or  to  those 
with  significant  loss.  For  nearly  lossless  dielectrics,  a  resonant  cavity  of  like  transverse 
geometry  as  the  three-dielectric  waveguide  was  used  to  investigate  the  physical 
phenomenon  and  verify  the  numerical  computations.  For  dielectrics  with  significant 
losses,  a  finite  waveguide  arrangement  through  a  shielded  enclosure  was  employed.  The 
results  of  both  measurement  techniques  compared  very  well  with  the  theoretical  predictions 
even  though  a  number  of  rough  approximations  had  to  be  made  to  recover  from  a  lack  of 
proper  equipment.  The  agreement  between  the  calculated  values  (the  resonant  frequencies 
of  the  cavity  and  the  attenuation  in  the  finite-guide  arrangement)  and  the  laboratory 
measurements  gives  credence  to  the  theoretical  and  numerical  work  in  this  thesis.  Coupling 
this  with  additional  verification  provided  by  other  published  work  for  the  two-dielectric 
waveguide  attests  to  the  accuracy  of  the  theoretical  development,  presented  next,  and  the 
subsequent  numerical  solutions  to  this  inhomogeneous  waveguide  arrangement. 
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CHAPTER  2 

THEORETICAL  FORMULATION 

The  primary  concern  of  this  chapter  is  the  theoretical  development  of  the  infinite 
inhomogeneous  waveguide.  The  particular  class  of  inhomogeneity  considered  here  is  a 
radial  variation  in  the  constitutive  parameters  forming  homogeneous,  lossy  annuli  with 
invariance  in  the  <J>  direction.  The  resulting  set  of  equations  is  too  complicated  to  solve 
analytically  since  the  expansion  functions  for  the  radial  variation  are  Bessel  functions; 
therefore,  numerical  techniques  are  employed.  The  principal  objective  of  this  numerical 
analysis  is  to  find  the  axial  wavenumber,  k2.  Once  this  has  been  determined,  all  other 
quantities  such  as  the  field  components,  attenuation,  and  power  flow  can  be  calculated. 
This  determination  of  is  based  on  an  eigensolution  concept  which  will  indicate  situations 
that  can  exist  and  not  necessarily  those  that  do  exist  in  a  physical  arrangement.  Further 
analysis  of  the  source  coupling  would  be  needed  for  this.  However,  since  a  sum  of 
solutions  is  also  a  solution  in  a  linear  treatment,  this  eigensolution  method  is  very  useful  in 
explaining  results  from  measurements  by  using  either  one  solution  or  perhaps  adding  a 
number  of  solutions. 

The  existence  of  hybrid  modes  in  this  type  of  inhomogeneous  guide  is  one  big 
factor  distinguishing  it  from  the  homogeneous  case.  A  hybrid  mode  is  one  that  has  both 
nonzero  axial  electric  and  magnetic  components.  This  is  in  contrast  to  the  homogeneous 
case  for  which  modes  of  vanishing  axiai  electric  or  magnetic  fields  can  be  constructed.  The 
method  of  mode  nomenclature  for  these  hybrid  modes  is  still  a  fuzzy  issue.  Many  authors 
have  put  forth  their  suggestions  but  no  method  has  won  sweeping  acceptance.  The  method 
chosen  here  is  somewhat  of  a  "hybrid"  in  that  Snitzer's  method  [18],  discussed  later,  is 
used  in  conjunction  with  a  mode-tracing  concept  outlined  by  Waldron  [19].  The  latter  topic 
is  discussed  fully  in  the  next  chapter.  It  is  stated  without  proof  by  Waldron  that  all  of  the 
modes  in  any  waveguide  have  a  1:1  correspondence  with  those  of  a  homogeneous,  pec 
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guide  of  appropriate  cross  section.  This  allows  a  method  of  nomenclature  based  on  the 
homogeneous  guide  as  a  limiting  case.  In  any  event,  each  mode  can  be  established  by  its 
axial  wavenumber.  Computation  of  this  wavenumber  is  possible  as  a  solution  to  a  set  of 
equations  constructed  from  first  principles.  This  formulation,  shown  next,  is  an  exact  one 
and  no  approximations  are  needed. 

2.1  Scalar-Wave  Function  Method 

The  cross  section  of  the  infinite  waveguide  having  multiple  concentric  circular 
dielectrics  is  shown  in  Figure  2.1.  (All  figures  and  tables  appear  at  the  end  of  each 
chapter.)  It  is  assumed  that  all  dielectrics  are  concentric  about  the  sheath  center  and  that  the 
constitutive  parameters  for  each  dielectric  are  constant.  Maxwell’s  equations  for  a  time- 
harmonic  field  in  a  homogeneous  space  are  given  by 


VxE  =  -joojiH  V  D  =  0 

VxH  =  jcoeE  V  B  =  0  . 


(2.1) 


The  time  dependence  eim  is  assumed  throughout  this  work  and  subsequently  suppressed. 
In  (2.1),  p.  is  the  permeability  of  the  medium  and  £  is  the  permittivity.  There  are  a  number 
of  ways  in  which  these  equations  may  be  combined  and  the  fields  solved  to  obtain  the 
characteristic  response  in  a  bounded  environment.  One  general  approach  is  to  stay  with  the 
six  field  components  and  write  solutions  in  terms  of  Ezand  Hz  [20].  Another  approach  is 
to  use  the  vector  potential  method  and  generate  solutions  to  a  two-dimensional  scalar 
problem  for  the  transverse  geometry  assuming  invariance  in  a  third  Cartesian  direction 
[21].  This  second  method  .  -hosen  here  for  two  major  reasons.  One,  it  is  elegant  in 
principle  since  all  the  field  co^nonents  for  a  mode  can  be  written  in  terms  of  one  scalar 


function  of  position.  Two,  some  interesting  properties  of  mode  designations  can  be  more 
clearly  delineated  when  scalar  wave  functions  are  used.  Whichever  method  is  chosen, 
identical  results  will  be  obtained.  The  only  differences  will  occur  in  the  intermediate  steps. 

What  follows  is  a  brief  review  of  the  vector  potential  method  to  solve  for  the 
characteristic  field  response  of  a  circular  cylindrical  waveguide.  This  formulation  is  needed 
to  introduce  the  methods  and  nomenclature  for  extending  the  normal,  homogeneous 
waveguide  methods  to  the  more  complicated  situation  with  q  concentric  annuli. 

Using  the  vector  potential  method  then,  these  equations  can  be  combined  and  the 
field  quantities  solved  yielding 


E  =  -VxF  -  jcopA  +  t^-V(V-A) 
J  **  joe 


H  =  VxA  -  jtoeF  +  t—~V (V'F) 

m 


(2.2) 


with  A  and  F  being  the  vector  magnetic  and  vector  electric  potentials,  respectively.  The 
vector  potentials  satisfy  the  vector  Helmholtz  equation 

V2T  +  k2T  =  0  (2.3) 

where  T  =  A  or  F  and  k2=to2p£.  Two  characteristic  field  responses,  namely.  Transverse 
Electric  (TE)  and  Transverse  Magnetic  (TM)  to  the  z  direction,  can  be  constructed  with 
appropriate  choices  of  these  vector  potentials.  Then  from  these  two  mode  types,  all  field 
patterns  in  the  guide  can  be  expressed  as  a  linear  superposition  of  them.  In  particular,  in 
the  i*  dielectric  let 


F,  =  0 


A,  =  zv^ 


(2.4) 
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where  vT  is  the  magnetic  scalar  wave  potential,  and  (2.2)  reduces  to 


E 

pi  jcoe^  dpdz 


i  a2< 

jCOCiP  d<pdz 


H„  = 


1^ 

p>~  p  dp 


(2.5) 


dp 


(  ^2 


j<*i 


fr^2K 


Vdz 


Hj.  =0 


with  being  the  intrinsic  wavenumber  of  the  medium.  Since  HZi=0,  the  field  is  TM  to  z. 
Likewise,  in  the  i*  dielectric  let 


Fj  =  zvf  A,  =  0 


where  V*  is  the  electric  scalar  wave  potential.  Then  (2.2)  reduces  to 


1  aV 

jcojij  5pdz 


P 


H, 


1 


d2Vf 


dpdz 


JWtiiP 


H 


1 

jam 


( zj. 


Uz2 


(2.6) 


Since  Ez,=0  in  this  case,  the  field  is  TE  to  z.  The  field  in  any  annulus  within  the  guide  can 
be  considered  as  a  linear  superposition  of  these  two  field  types. 
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Also,  since  the  Vi's  are  chosen  as  Cartesian  components  of  the  vector  potentials, 
the  Vi's  necessarily  satisfy  the  scalar  Helmholtz  equation 


V2Vi  +  k,2Vi  =  0  .  (2.7) 

In  addition,  the  geometry  inside  the  guide  remains  separable  with  the  inclusion  of  these 
concentric  dielectrics.  Therefore,  it  is  quite  easy  to  construct  solutions  to  the  scalar 
Helmholtz  equation  in  each  dielectric.  Solutions  to  this  homogeneous  partial  differential 
equation  can  be  found  using  the  separation  of  variables  method.  Accordingly  then, 

Vi  =  P(p)<D(<|>)Z(z)  .  (2.8) 

By  substituting  this  equation  into  (2.7)  and  expressing  the  Laplacian  in  cylindrical 
coordinates,  (2.7)  can  be  separated  into  three  equations,  each  a  scalar  homogeneous 
ordinary  differential  equation  in  one  cylindrical  coordinate  only.  For  the  circular  cylindrical 
geometry  here,  the  general  form  of  v,  will  be 

Vi  =  BnOcp.pMntlOhOe^)  (2.9) 

where  Bn  is  some  linear  combination  of  solutions  to  Bessel's  equation,  h(n<f>)  is  some  linear 
combination  of  the  trigonometric  functions  sin(n<f>)  and  cos(n(j>),  and  hfk^z)  represents  the 
longitudinal  variation  of  the  propagating  wave.  For  a  wave  propagating  in  this  axially 
invariant  environment 


h(kzz)  =  e-^2 


(2.10) 
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is  a  correct  choice.  Here,  lq  represents  the  longitudinal  wavenumber.  In  a  transversely 
bounded  environment,  the  Bn  function  will  be  a  linear  combination  of  two  independent 
solutions  to  Bessel's  equation  of  the  form 


Bn(kpp)  =  [aiJn(kPip)+biNn(kp.p)l  (2.11) 

with  Jn  and  Nn  the  Bessel  functions  of  the  first  and  second  kinds,  respectively,  and  ai  and 
bL  being  constants.  For  any  source-free  region  containing  the  origin,  b  will  be  zero.  For  a 
physical  system,  one  would  expect  the  fields  to  be  single  valued.  This  implies  that  n 
should  be  chosen  as  an  integer.  Finally,  for  the  to  be  a  valid  solution  to  (2.7)  the 
separation  relation 

k2  =  co2^  =  k2  +  k^  (2.12) 

must  be  satisfied.  By  phase  matching  considerations  of  the  tangential  fields,  1^  will  be  the 
same  throughout  the  cross  section  of  the  guide. 

2.2  Multiple  Concentric  Dielectrics 

To  formulate  the  problem  of  q  concentric  dielectrics,  the  Vj/j's  are  constructed  in 
view  of  (2.9)  and  the  ensuing  discussion.  In  each  dielectric  then. 


vT  =  CnJn(kPip)cos(n0)e~jklZ 
Vi  =  CnJ^kp^lsinlnple*^ 


(2.13) 


yg  =  [C2iJ„(kp2p)  +  C22Nn(kPip)lcos(n«)e* 
V2  =  iCaaJ.aw*) +  C24Nn(kp2p)]sin(n4))e_jklZ 


VT  =  fCtiJn(kptp)  +  Ct2Nn(kpp)]cos(n«{>)e_jk*z 
Vf  =  [CaJ^kpp)  +  Ct4Nn(kpp)]sin(n<t>)e--ikl2 

<  =  [CqlJn(kpqp)  +  Cq2Nn(kPqp)]cos(n<{>)e-^z 

<  =  [Cq3Jn(kpqP)  +  Cq4Nn(kPqp)]sin(n<|))e"j1^  . 

The  coefficients  C„  are  unknowns  at  this  stage  and  will  be  solved  for  in  the  next  chapter. 
The  index  r  refers  to  the  layer  and  s  to  the  coefficient  of  the  specific  radial  variation  in  the 
expansion  (an  integer  from  1  to  4).  C12  and  C14  are  both  chosen  equal  to  zero  since  the 
Neumann  function,  Nn,  is  infinite  at  the  origin  and  there  is  no  reason  to  expect  singularities 
in  the  field  there.  By  the  azimuthal  symmetry  of  this  waveguide  arrangement,  there  is  a 
degeneracy  in  the  <}>  variation  except  for  the  modes  having  n=0  [7].  For  simplicity,  only 
one  of  these  is  shown  here. 

The  boundary  conditions  for  this  problem  (that  the  tangential  fields  are  continuous 
across  each  dielectric  boundary)  must  be  applied  in  order  to  evaluate  the  coefficients.  At 
each  dielectric  interface  there  will  be  four  equations  of  continuity —  one  equation  for  each 
of  Ep,  Ez,  Hp  and  Hz.  This  gives  a  total  of  (4q-4)  equations.  The  total  number  of 

unknowns  is  (4q-2)  due  to  the  fact  that  C12  and  C14  have  been  evaluated  on  physical 
grounds.  Two  more  equations  are  needed.  These,  of  course,  are  provided  by  the  guiding 
structure  which  determines  the  major  propagating  characteristics  of  the  system.  Assuming 
in  this  case  a  pec,  then  Ep=0  ana  Ez=0  at  p=b.  Therefore,  the  total  number  of  equations 

and  unknowns  is  equal  at  (4q-2). 
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As  mentioned  previously,  the  total  field  in  each  annulus  will  be  a  linear 
superposition  of  TE  and  TM  fields  to  z.  To  apply  the  boundary  conditions  then,  the  fields 
must  be  computed  according  to  Equations  (2.5)  and  (2.6).  Proceeding,  let 

VT  =  V^hm(n4>)  =  {[C.jVkpp)  +  Ci2Nn(kPip)}e_J^}cos(n<})) 

and  (2.14) 

=  vfhe(n<J>)  =  {[Ci3Jn(kp.p)  +  Ci4Nn(kp  p)]e~jklZ}sin(n<J))  . 


Using  this  notation. 


*p. 

e^  =  — vf1 

*1  JCOfij Yl 


TT  **  C 

H  -  - - w. 


(2.15) 


(2.16) 


and 


jco^p  ~W  "3p” 

E*. =  +  kPiV^  h<:(n<i>)  • 


Also, 


H 

H 


-jkz  fyj  fyT 

jCD^p  dp  dp 


(2.18) 
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The  above  primed  quantities  denote  differentiation  with  respect  to  p  as 

vl" '  =  [C„  V(kPip)  +  Cj2Ntl,(kp.p)le^k*z 

V? '  =  [Cl3I„'(kpP)  +  C^N/d^pIle-^1  . 


Equating  tangential  components  of  E  and  H  at  p=a  gives 


(i)  for  Ez  — 


at  p=a 


then. 


kP(hD 

jcae<t_i) 


or, 


+  Qt-l)2^n(kp(tl)a)] 

~  ^p,^(i-l)tQl^n^pta^  +  Q2^n^pta^  =  0 


(ii)  for  Hz  — 

Hz<I_1)  =  Hz,  ^  P  =  a 

then, 

Wo  ^  *e 

j%1)'|,y 

or, 


+  Qt-D4Nn(kp(il)a)] 

—  ^),M,(t-i)fCt3Jn(kpta)  +  Ct4Nn(kp  a)]  =  0 


(2.19) 


(2.20) 


(2.21) 


(iii)  for  E0  - 

then, 

or. 


(iv)  for  — 

then. 


at  p  =  a 


coe(t-i)a 


Va-i)he(n<{>)  +  kP([_|)Va-i)he(n<j>) 


nkz 

coe^ 


vThe(n<}>)  +  he(n<j>) 


=  0 


coe^a  ^C(‘-i)iJn^p(l.,)a)  +  Qt-mNn^Sv,^ 

+  ^P(t_t)^h-l)3^n  ^(t.,)a^  +  Qt-l)4^n  (kP(,.0a)] 


nkj 

toe^a 


[CtlJ„(kpa)  +  CgN^kp  a)] 


-  kp/CoVfkp  a)  +  Ct4Nn'(kp  a)]  =  0 


n^a-D  H«t 


at  p  =  a 


^  ¥a-i)hm(n4>)  +  kP([_i)y(tnl1')hm(n0) 


(04<t_i)a 
nkj 


cot^a 


\|/fhm(n0)  +  kp  ytm  hm(n<J>) 


=  0 


^Pa-J^t-Di^n  (kP(tI)a)  +  Qc-iteNn'Ckp^a)] 


nk* 


+  a)|i(t_[)a  [C(t-i»J"(kP(t-oa)  +  C(t-D4Nn(kPa.0a)] 


or. 
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-  kpfcrtVfcp  a)  +  CaN.Okp  a)J 


co^ta 


[Ct3Jn(kpa)  +  Ct4Nn(kpa)]=0  . 


The  remaining  boundary  condition  to  apply  is  that  for  the  guiding  structure  itself.  For  a 
pec,  the  tangential  electric  fields  Ep  and  Ez  must  equal  zero  at  p=b.  This  implies  from 

considering  Equations  (2.13),  (2.15)  and  (2.17) 


<V„«S,W  +  CqjNjjOCp^b)  =  0 

and  (2.24) 

C^V(kpW  +  b)  =  0  . 


Equations  (2.20)  through  (2.24)  form  the  primary  means  by  which  this  inhomogeneous 
waveguide  problem  will  be  solved.  Collectively,  these  equations  may  be  considered  to 
form  a  (4q-2)  matrix  on  the  "variables"  C„.  Since  there  are  no  forcing  terms,  these 
equations  can  be  written  in  the  matrix  form 

[Al[C]  -  0  (2.25) 

where  [A]  is  a  (4q-2)  square  matrix  and  (C]  is  a  (4q-2)  column  matrix.  For  the  nontrivial 
cases,  (2.25)  is  valid  only  for 


det[A]  =  0  .  (2.26) 

The  variable  that  will  be  solved  for  in  order  to  enforce  (2.26)  will,  of  course,  be  kz.  Once 
k,,  is  known,  all  other  variables  and  fields  can  be  computed  as  will  be  shown  later. 
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2.3  Two-Dielectric  Matrix  Construction 

It  may  be  instructive  at  this  point  to  actually  construct  the  matrix  [A]  for  the  case 
having  two  dielectrics.  This  matrix  will  then  be  compared  to  other  work  as  a  verification. 
Let  the  inner  dielectric  have  constitutive  parameters  jj.],  £i  with  radius  a,  and  the  outer  p.2, 
e2  with  radius  b.  The  guiding  structure  will  be  assumed  a  pec.  The  matrix  [A]  is  then 
constructed  using  Equations  (2.20)  through  (2.24)  giving 


kJ^E! 

0 

kp^Ej 

^iE4  0 

0 

0 

kpj|42Ei 

0 

0  kp^lE2 

kp2^iE4 

— E 
coe^1 

kpEi' 

nkz  _ 
(ue2a  2 

nlc, 

— ^-E4  k, 
oje2a  4  1 

kpE4' 

kp.Et' 

nk, 

——•Et 
top^a  1 

kp2E2' 

nk, 

kD  E4'  — 2-E2 
P2  4  cop.2a  1 

nk, 

cop2a  4 

0 

0 

e2 

E4  0 

0 

0 

0 

0 

0  E2' 

e4' 

(2.27) 


where 


and 


Ei=J„(kpa) 
E2  =  J„(kp2a) 


£3  =  ^(1^) 
E4  =  Nn(kp2a) 


[C]  = 


C13 

C21 

C22 

C23 

.C24_ 


with 


18 


(A][C]  =  0  . 

This  matrix  is  in  agreement  with  Harrington  [21].  Here,  the  matrix  could  be  made  smaller 
in  size,  from  (4q-2)  down  to  (4q-4)  through  elimination  of  the  last  two  rows  by  combining 
the  third  and  fourth  columns  and  the  fifth  and  sixth  columns.  However,  to  be  more  general 
and  allow  for  easier  computer  programming  with  the  inclusion  of  other  guiding  structure 
boundary  conditions,  this  form  is  preferred. 

2.4  Hybrid  Modes 

Examining  matrix  [A],  when  n=0,  coefficients  Cn,  and  are  not  related  to 
coefficients  C13,  Q3  and  C24  and  the  field  separates  into  modes  TE  and  TM  to  z  [21], 
From  (2.14)  then,  the  y*  's  are  not  related,  or  coupled,  to  the  y™  's.  This  is  precisely 
what  is  meant  by  TE  or  TM  "modes."  That  is,  there  exists  a  relationship  between  five  of 
the  field  components  describing  an  allowed  field  pattern  with  the  sixth  field  component, 
either  Ej,  or  Hj,  equal  to  zero.  But  there  is  no  relationship  between  the  coefficients  and 
hence  the  fields  of  the  two  mode  types.  In  all  other  cases,  where  n*0  and  1^*0,  the  fields 
are  neither  TE  nor  TM,  but  are  hybrid  modes.  This  means  that  each  mode  in  the  guide  will 
have,  in  general,  both  nonzero  Ez  and  H^. 

The  main  emphasis  here  is  that  the  field  pattern  in  each  dielectric  is  considered  to  be 
a  superposition  of  two  "basis  fields":  the  TE  and  TM  modes.  This  is  done  in  each  dielectric 
and  the  \y; 's  are  constructed  as  per  separation  of  variables  of  the  scalar  Helmholtz 
equation.  The  boundary  conditions  are  enforced  on  the  system  and  the  expansion 
coefficients  are  solved  for.  If  the  expansion  coefficients  of  one  mode  type  are  not  related 
to,  or  not  coupled  to,  the  expansion  coefficients  of  the  other  mode  type,  then  TE  and  TM 
fields,  or  modes,  can  exist  in  the  guide.  Herein  lies  the  main  point.  This  means  that  the 
field  can  have  its  pattern  described  by  just  the  TE  field  relationships  or  just  the  TM  field 
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relationships.  However,  if  the  guide  is  inhomogeneous,  the  coefficients  of  the  TE  and  TM 
modes  become  interdependent  Then  the  field  must  be  written  as  a  superposition  of  TE  and 
TM  field  patterns  and  cannot  be  separated.  The  expansion  coefficients  of  the  two  mode 
types  become  coupled  in  this  case. 

Another  way  to  look  at  this  situation  and  investigate  this  coupling  of  coefficients 
concept  is  to  go  back  and  consider  Maxwell's  equations  in  (2.1).  Not  every  solution  to 
(2.3)  (with  T  =  E  or  H)  will  satisfy  Maxwell's  equations.  One  also  needs  to  enforce  the 
divergence  relations  in  (2.1).  Expanding  these  out  in  cylindrical  coordinates  and 
substituting  the  linear  superposition  of  (2.5)  and  (2.6)  for  the  field  components  in  a  source- 
free  homogeneous  medium  give  in  the  i*  dielectric 


V-CjEj  =  0  => 


i  a 

favrl 

\  dvf 

i  a 

i 

I'sV) 

1  avf 

p  dp 

— 
.  V 

l  apaz ) 

I’6-!*-. 

.m 

l  apaz  J 

1  +  e‘‘3p". 

+ 


a 


rk2 

V  JO)  Yl 


=  0 


(2.28) 


i  a 

p  f 

aV| 

i  a 

i  | 

f dVr  T 

a 

f  J 

P  dp 

.  apaz  / 

.  +  P  dp 

.]wp 

l  apaz  J 

.  +  dz 

JO)  Yl  ) 

(2.29) 


which  is  only  a  function  of  Likewise, 
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which  is  only  a  function  of  y*.  What  has  tacit’ y  been  assumed  here  is  that 


1  b  ayf  £j  a V  }_b_  3vT  m aV* 

p  dp  dp  p  dpdp  °r  pdp^1  dp  p  dpdp 


This  is  true  only  if  these  partial  derivatives  are  continuous  and  the  medium  homogeneous. 
All  of  the  functions  in  (2.32)  are  discontinuous  at  p=0.  Otherwise,  they  are  continuous  in 
a  homogeneous  space  so  that  Equations  (2.29)  and  (2.31)  are  valid,  and  the  coefficient  sets 
(Qi.Cjj)  can  be  determined  independently  of  the  sets  (C^.C,^.  Hence,  TE/TM  modes  are 
possible. 

However,  in  an  inhomogeneous  guide  this  is  no  longer  the  case.  Consider  from 
(2.28)  and  (2.30) 


i  a  ¥  ti  dV'  ,  i  a  avT  Hi  aS'T 

“d  pIp^TT  *  <233) 

since  in  general,  £t=f(p)  and  p.i=g(p).  (No  variation  of  Ejand  in  <j>  is  assumed  in  this 
study.  If  this  is  not  the  case,  modal  separation  of  fields  TE  and  TM  to  z  would  not  be 
possible  and  only  hybrid  modes  would  exit.)  In  order  for  the  interchanging  of  the  partials 
with  respect  to  p  and  <j>  to  be  a  valid  mathematical  operation,  the  mixed  partials  of  \j/;  must 
be  continuous.  In  an  inhomogeneous  guide  this  will  generally  not  be  true.  The  y;'s  will 
be  discontinuous.  Consider  Figures  2.2  and  2.3  which  are  graphs  of  the  Vi’s  for  a  2.005 
in  radius  pec  waveguide  having  a  0.1908  in  thick  Teflon  tube  (£,.=2.1  F/m)  of  inner  radius 
0.5008  in  filled  with  lossy  methanol  (£,=23.9  -  j  15.296  F/m)  at  3  GHz.  (The  permittivity 


21 


values  are  from  von  Hippel  [22],  and  the  xj/j's  and  their  derivatives  will  be  solved  for  in 
Chapters  3  and  4.)  In  both  cases,  the  x^'s  and  their  derivatives  are  discontinuous,  within 
the  discretization  of  the  plot,  across  each  material  boundary  but  continuous  within  each 
annulus.  This  is  in  accordance  with  the  previous  discussion.  These  xj/j's  are  discontinuous 
functions  in  p  and  interchanging  partial  derivatives  is  not  a  valid  mathematical  operation. 
Therefore,  the  two  coefficient  sets  are  related  as  in  (2.28)  and  (2.30).  The  manner  in 
which  the  coefficients  of  the  two  mode  types  are  coupled  is  not  as  important  as  the  fact  that 
they  are  just  related.  From  this  information  alone  it  is  sufficient  to  deduce  that  TE  and  TM 
fields  cannot  exist  in  the  guide,  but  field  modes  having  all  six  components  each  must  exist, 
i.e.,  hybrid  modes.  Also  evident  in  these  graphs  is  at  p=b 


This  is  consistent  with  the  boundary  conditions  imposed  on  the  system  for  and  EZ3. 
Evidenced  in  this  example  is  perhaps  the  biggest  difference  between  the  scalar  wave 
potential  and  field  approaches  to  solving  these  EM  problems.  In  large,  the  xy's  will  be 
discontinuous  functions,  whereas  the  fields  will  be  continuous  except  for  bounded-charge 
induced  discontinuities. 

Another  interesting  case  is  that  of  1^=0.  For  ail  n,  the  hybrid  mode  separation  into 
modes  TE  and  TM  also  form  identically  to  the  cases  when  n=0.  As  an  example,  the  matrix 
(2.27)  constructed  for  two  dielectrics  is  informative.  Here,  n  and  kj  occur  in  pairs.  When 
kz=0  the  same  behavior  as  for  n=0  will  occur,  that  is,  hybrid  modal  separation  into  TE  and 
TM  modes.  However,  this  separation  is  for  all  modes  at  cutoff  and  not  just  the  n=0 
modes.  This  gives  rise  to  another  modal  designation  scheme,  namely,  as  the  mode  tends  to 
cutoff,  the  correspondence  between  the  homogeneous  cases  can  be  correlated.  Although  it 
has  been  reported  that  this  method  of  classification  agrees  with  the  limiting  case  as  the 
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guide  becomes  homogeneous  [23],  this  is  not  always  true.  If  no  backward-wave  region 
(Section  3.3)  is  present  in  the  dispersion  graph,  the  above  statement  is  correct.  As  will  be 
seen  in  Sections  (3.5)  and  (4.2),  with  a  backward-wave  region  present,  a  small  amount  of 
loss  must  be  introduced  into  the  dielectrics  to  choose  the  appropriate  eigenvalue  and 
corresponding  mode  designation.  It  was  found  from  the  results  after  applying  this 
technique  that  these  two  schemes  do  not  always  agree  if  cutoff  is  defined  by  the  frequency 
where  kj  is  zero.  Garricoats  and  Taylor  [14]  note  something  similar  to  this  but  place  more 
restrictions  on  the  rod  insert  ^  value  which  yields  this  effect  than  is  really  necessary. 

For  more  than  two  dielectrics,  when  either  n  or  kj  is  zero,  the  modes  become  TE 
and  TM  in  the  same  fashion  as  with  two  dielectrics.  As  each  new  dielectric  is  added,  the 
original  coefficient  sets  of  the  two-dielectric  matrix  remain  uncoupled  and  four  new 
columns  are  added  in  which  the  new  coefficient  sets  themselves  arc  uncoupled.  These 
columns  are  added  in  virtue  of  Equations  (2.20)  through  (2.23).  From  (2.20)  and  (2.21)  it 
is  seen  that  the  z  components  of  the  fields  do  not  provide  this  coupling.  It  is  only  the  <J> 
components  which  contribute.  With  an  eye  on  (2.22)  and  (2.23),  when  n  or  1^  =0  the 
(Cti,Ct2)  become  uncoupled  from  all  other  coefficients,  as  do  (C^.C^).  The  Ct 
coefficients  have  the  same  form  as  the  's  so  the  two  sets  of  coefficients  in  the 
previous  four  columns  remain  uncoupled  with  the  addition  of  four  new  columns. 
Similarly,  when  n  and  kz  are  not  zero,  the  addition  of  four  new  columns  for  each  new 
dielectric  interface  will  not  have  any  decoupling  effect  or  leave  any  other  marked 
impression  on  the  matrix  [A]  other  than  increasing  the  size.  The  work  on  modal 
designations  that  has  been  applied  to  two  dielectrics  here  and  in  the  literature  can  be 
extended  to  three  or  more  dielectrics  without  modification. 
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2.5  Modal  Nomenclature  and  Designation 

Quite  surprisingly,  there  is  no  universal  scheme  of  nomenclature  for  modes  in 
dielectrically-loaded  waveguides.  There  have  been  many  attempts  and  Bruno  and  Bridges 
[6]  give  an  excellent  history  of  the  major  ones  including  a  new  suggestion  of  their  own. 
Beam  had  one  of  the  earliest  methods,  which  was  based  on  the  relative  contributions  of  Ez 
and  Hz  to  a  transverse  field  component  [24].  If  Hj  made  a  larger  contribution,  the  mode 
was  designated  HE,  and  if  Ez  larger,  EH.  The  fundamental  was  assigned  EHn.  This 
scheme  was  reported  to  be  arbitrary  since  it  depends  on  which  transverse  component  was 
observed  and  how  far  the  mode  is  from  cutoff.  Snitzer  devised  a  new  scheme  by  which 
most  modem  day  methods  are  based  [18].  Although  Snitzer’s  work  was  with  a  dielectric- 
rod  waveguide,  the  extension  to  the  bounded  guide  is  trivial.  This  method  considers  the 
HEn  mode  fundamental  by  both  definition  and  common  usage  at  the  time  (1961).  All 
modes  that  have  the  same  sign  as  the  fundamental  for  the  ratio  of  the  coefficients  of  EZ[  and 
HZl,  namely,  C^/C^,  are  designated  HE.  All  other  modes  with  different  signs  are 
designated  EH.  Note  that  the  use  of  HE  and  EH  here  is  just  opposite  to  that  of  Beam's. 
From  numerical  work  in  Section  4.2,  this  scheme  works  well  for  modes  above  cutoff  but 
apparently  does  not  apply  to  evanescent  modes  since  no  consistent  relationship  among  the 
coefficients  could  be  found. 

An  efficient  and  unambiguous  method  of  classification,  which  can  be  used  in 
conjunction  with  Snitzer's  scheme,  is  to  start  with  the  homogeneous  case,  either  a  TE  or 
TM  mode,  and  gradually  increase  (or  decrease)  the  radii  of  the  dielectrics  within  the  guide 
and  trace  kz.  Thereby,  the  modes  HE  and  EH  can  be  associated  with  TE  and  TM  modes, 
respectively,  of  the  homogeneous  waveguide.  This  is  the  Correspondence  Idea  of  Waldron 
[19]  from  which  a  1:1  correspondence  of  all  modes  in  the  dielectric- filled  guide  is 
maintained  with  those  of  the  homogeneous  pec  guide  as  all  nonideal  characteristics  of  the 


former  approach  those  of  the  latter.  For  example,  the  dielectric  annuli  radii  approach  zero, 
the  eri  approach  one,  or  the  losses  in  the  dielectrics  or  waveguide  walls  approach  zero. 
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2.6  Nonzero  Wall  Losses 

If  in  a  practical  situation,  the  metal  losses  were  significant  compared  to  the  dielectric 
losses,  the  formulation  for  this  situation  can  be  accurately  approximated  by  considering  the 
metal  to  be  infinite  in  extent  [25].  If  the  skin  depth  is  small  with  respect  to  the  guide  wall 
thickness,  this  is  a  good  approximation.  Then  instead  of  Jn  and  Nn  functions,  I„  and  K„ 
can  be  used  where  In  is  the  modified  Bessel  function  of  the  first  kind  and  is  the 
modified  Bessel  function  of  the  second  kind.  Since  exhibits  e_kp  asymptotic  behavior 
while  ^  has  e1^  behavior,  the  coefficients  multiplying  the  Ij,  functions  must  be  set  equal  to 
zero.  That  is 

Vq*  =  [Cql^p^P)  +  Cq2Kn(kpqp)]c  ^  (2.34) 

¥q  =  [CqaU^P)  +  Cq4Kn(kPqp)]e-^z 

where  Cql=Cq3=0  and  there  are  (q-1)  dielectrics  in  the  waveguide.  No  other  boundary 
conditions  need  be  applied  to  this  system  and  the  size  of  the  [A]  matrix  is  reduced  to 
(4q-4).  This  example  shows  that  the  analysis  used  in  this  chapter  is  general  enough  that  it 
is  relatively  easy  to  modify  the  equations  and  corresponding  matrices  to  account  for  other 
waveguiding  boundaries,  not  just  pec  waveguide  structures. 

To  actually  find  solutions  to  these  radially  inhomogeneous  waveguide  problems, 
the  characteristic  equation  (2.26)  is  solved  yielding  the  eigenvalues  This  matrix  [A],  as 
noted  earlier,  is  composed  of  higher  transcendental  functions  and  in  general  cannot  be 
solved  analytically.  Instead,  numerical  techniques  are  employed,  which  is  the  subject  of 


25 


the  next  chapter.  It  should  be  noted  in  passing  that  certain  special  cases  may  be  solved 
analytically  (for  example,  n=  0  or  1  and  two  dielectrics)  and  using  tables  of  Bessel 
functions  to  obtain  numerical  results.  However,  if  lossy  materials  are  considered  or  higher 
modes  are  needed,  solutions  are  not  possible  using  only  analytical  techniques  and  tables  of 
Bessel  functions,  even  for  these  special  cases. 


Magnetic  scalar  potential, and  its  derivative 
Teflon  pipe  filled  with  methanol  at  3  GHz. 
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CHAPTER  3 

NUMERICAL  IMPLEMENTATION 

The  characteristic  equations  constructed  in  Chapter  2  are  too  complicated  to  solve 
analytically,  and  require  a  numerical  treatment  instead.  This  chapter  is  concerned  with 
computation  of  these  solutions  through  implementation  of  a  root-finding  routine  to  search 
for  the  roots  of  these  characteristic  equations.  A  group  of  five  Fortran  computer  programs 
was  written  to  accomplish  this  task  plus  many  others,  including  the  calculation  of  k^  versus 
0)  plots,  kz  versus  relative  radius  plots,  transverse  power  flow  and  more.  With  space 
limitations  present,  only  one  code,  MDCW,  is  listed  here.  The  other  codes  are  basically 
derivatives  of  this  main  one. 

The  use  of  these  programs  and  interpretation  of  the  ensuing  data  will  be 
demonstrated  first  for  the  two-dielectric  case.  This  dielectric-rod  insert  guide  case  has  been 
thoroughly  investigated  in  the  past  as  was  noted  in  Chapter  1.  Comparison  of  the  results 
here  with  those  in  earlier  published  work  will  lend  some  validity  to  the  theoretical 
development  and  confirm  the  accuracy  of  the  computer  programs.  In  addition  to  this,  some 
extensions  of  this  earlier  work  will  be  presented.  Most  notable  of  these  is  in  the  area  of 
power  flow,  where  it  will  be  shown  that  the  idea  of  a  backward  wave  becomes  more 
sharply  defined  for  dielectrically-loaded  waveguides. 

Much  emphasis  in  this  chapter  is  placed  on  the  two  hybrid  modes  HEt  j  and  EHn. 
The  reasons  for  this  are  that  the  HEn  mode  is  always  the  fundamental  mode  (has  the 
lowest  cutoff  frequency)  and  that  the  EHn  mode  was  found  to  exist  in  some  finite 
waveguide  experiments  conducted  in  the  laboratory.  Some  astonishing  modal  conversions 
between  these  two  modes  were  observed  in  these  experiments.  Hence,  for  a  comparison  of 
the  theoretical  predictions  and  the  experimentally  measured  results,  these  modes  are 
emphasized. 
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3.1  Root  Finding  by  Muller's  Method 

Solutions  to  (2.26)  must,  in  general,  be  obtained  numerically.  This  equation  can  be 
thought  of  as  a  one-dimensional  function  in  the  argument 

fOO  =  0  .  (3.1) 

However,  if  losses  are  present  or  in  some  fashion  kj  becomes  complex  such  as 

kz  =  p-ja,  (3.2) 

then  the  one  dimensionality  of  the  function  is  lost  and  solutions  to  (2.26)  must  be  found  in 
the  complex  plane.  Probably  the  most  popular  numerical  method  for  solving  a  generic 
problem  such  as  (3.1)  is  the  Newton- Raphson  method.  The  difficulty  in  applying  this 
method,  however,  is  that  the  derivative  of  the  function  must  be  computed.  Since  the 
derivative  is  unknown  analytically  and  numerical  differentiation  is  not  a  stable  operation, 
the  Newton-Raphson  method  will  not  be  used.  Instead,  the  roots  of  (3.1)  will  be  found 
using  Muller's  method.  This  method  is  a  generalization  of  the  widely  known  secant 
method,  but  an  inverse  quadratic  interpolation  is  used  instead  of  a  two-point  linear 
interpolation  [26].  Muller's  method  can  be  derived  by  considering  the  quadratic  function 
[27] 


f(x)  =  a(x-x2)  +  b(x-x2)  +  c 


(3.3) 


The  coefficients  a,  b,  and  c  can  be  evaluated  with  the  specification  of  three  initial  guesses, 
x0,  xb  and  x2,  for  the  root.  The  computed  root  of  this  quadratic  function  is 
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x3  =  x2  -  (x2-xj) 


2c 

b±V  b2-4ac 


(3.4) 


with 

x2-xi 

q  - - 

xi-xo 

a  3  qf(x2)  -  qd+qMxi)  +  q2f(xo) 
b  s  (2q+l)f(x2)  -  (l+q)2f(x!)  +  q2f(x0) 
c  s  (l+q)f(x2) 

A  visual  image  of  this  quadratic  fit  and  the  implementation  of  this  method  are  shown  in 
Figure  3.1.  The  zeros  of  this  unique  quadratic  are  calculated  and  the  next  point  for  another 
quadratic  fit  is  located  using  xlf  x2,  and  now,  x3.  This  process  is  repeated  until  the  zero  of 
the  exact  function  is  found  within  a  specified  tolerance.  The  sign  in  (3.4)  is  chosen  to 
make  the  magnitude  of  the  denominator  the  largest.  The  two  main  advantages  of  this 
method  are  that,  one,  complex  roots  may  be  found  due  to  the  quadratic  nature  of  the 
interpolation  and,  two,  the  root  does  not  necessarily  have  to  be  bracketed  with  the  initial 
guesses  in  order  to  converge.  In  fact,  just  as  with  the  secant  method,  the  root  may  not 
remain  bracketed  after  the  searching  process  begins  even  if  it  initially  was.  This  advantage 
may  also  be  a  disadvantage,  since  for  certain  functions,  Muller’s  method  may  actually 
diverge  from  the  root  depending  upon  the  behavior  of  the  exact  function  near  the  root.  All- 
in-all,  Muller's  method  is  very  useful  and  accurate  for  most  applications  of  this  type.  It  has 
its  own  shortcomings  in  certain  instances,  most  of  which  can  be  rectified  by  simply 
restarting  the  search  routine  and  changing  the  initial  guesses. 
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3.2  MDCW  Computer  Program 

Muller's  method  forms  the  basis  of  the  computer  program  "MDCW"  (Multiple 
Dielectric  Circular  Waveguide)  which  is  used,  among  other  things,  to  solve  the 
characteristic  equation  (2.26).  (A  listing  of  the  Fortran  source  code  is  given  in  the 
Appendix  along  with  a  brief  description  of  the  program.  This  program  is  the  basic 
computer  code  used  in  the  study  of  this  waveguide  problem.  Most  of  the  other  programs 
written  and  used  in  this  work  are  derivatives  in  whole  or  in  part  of  the  philosophy  and 
methods  employed  here.)  In  order  to  search  for  the  roots  of  (2.26),  the  determinant  of  the 
matrix  must  be  computed  repeatedly.  That  is,  in  each  iteration  of  this  root-finding  scheme, 
the  function,  f(k2),  being  evaluated  is  the  determinant  of  [A].  An  efficient  method  for 
finding  the  determinant  of  a  matrix  is  by  LU  decomposition  [26].  This  is  a  tridiagonalizing 
method  whereby  the  determinant  is  computed  as  the  product  of  the  diagonal  elements  in 
either  an  upper  or  lower  triangular  matrix.  The  LU  decomposition  can  also  be  used  to 
solve  a  set  of  linear  equations  by  back  substitution  which  will  be  of  value  later  on  when 
computing  the  coefficients  C„.  This  method  then  is  efficient  in  that  it  can  perform  double 
duty  so  to  speak. 

The  iteration  process  searching  for  the  root  kz  that  satisfies  the  characteristic 
equation  is  terminated  upon  satisfaction  of  some  predetermined  convergence  criterion.  The 
criterion  chosen  here  is  for  two  kz  values  to  differ  by  less  than  some  user  specified 
accuracy.  Typically  10*8  or  10'9  is  used.  Although  these  are  quite  small  numbers, 
Muller's  method  converges  extremely  fast  for  the  functions  encountered  here.  In  addition, 
specifying  small  tolerances  such  as  these  can  actually  speed  up  convergence  in  some 
instances  as  when  this  program  is  automated  for  kz  versus  co  plots.  Calling  these 
"tolerances"  an  "accuracy”  can  be  somewhat  misleading  since  convergence  to  incorrect 
answers  is  also  possible.  Once  a  certain  root  has  been  located,  other  factors  must  be  taken 
into  account  to  verify  that  this  is  a  valid  root  and  the  one  being  sought.  These  factors 


include  the  value  of  the  determinant,  the  values  of  the  coefficients  and  kpj's,  and  the 
resulting  field  patterns  which  can  be  computed  once  the  coefficients  have  been  solved. 
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After  the  kj  value  has  been  found,  the  corresponding  [C]  vector  (cf.  (2. 251)  can  be 
calculated.  Since  det[A]=0  is  the  root-finding  criterion,  the  rows  and  columns  of  [A]  are 
not  all  linearly  independent.  In  general,  with  n=(4q-2),  the  rank  of  [A]  is  (n-1).  Since 
there  are  n  unknowns  and  now  (n-1)  linearly  independent  equations,  one  component  of  the 
[C]  vector  can  be  chosen  arbitrarily.  In  this  case,  Cn  will  be  chosen  equal  to  one.  (The 
reason  Cn  is  chosen  and  not  another  component  is  important  and  will  become  evident 
when  power  flow  is  discussed.)  An  inhomogeneous  set  of  equations  results  of  the  form 


a22  ' 

a2n 

^13 

a21 

.^2 

ann_ 
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where  the  a's  are  elements  of  the  [A]  matrix.  This  system  of  equations  is  solved  giving  the 
Qs's  using  LU  decomposition  plus  back  substitution  [26]. 

The  final  step  in  computing  the  fields  from  (2.5)  and  (2.6)  is  to  substitute  the  C^'s 
into  (2.13).  This  then  will  give  all  of  the  fields  throughout  the  guide.  As  has  been 
mentioned  in  Section  2.4,  the  fields  within  the  guide  will  be  a  linear  superposition  of  fields 
TE  and  TM  to  z. 

One  principal  requirement  for  this  program  was  that  it  be  able  to  run  on  a  personal 
computer  (PC).  Using  Microsoft's  Fortran  v4.01  and  double  precision  real  and  complex 
numbers  (64  bits  of  precision  for  each  of  the  real  and  imaginary  parts),  one  total  iteration 
for  three  dielectrics  takes  approximately  2  sec  on  a  12  MHz  clocked  286  AT  compatible  PC 
with  an  8  MHz  80287  math  coprocessor.  On  a  386  Compaq  running  at  16  MHz  with  a  16 
MHz  80387  math  coprocessor,  the  iteration  time  is  reduced  to  approximately  1  sec. 
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Typically,  for  initial  guesses  within  5  to  10%  of  the  final  answer,  4  to  8  iterations  are 
usually  needed  to  converge.  The  total  times  needed  to  find  a  solution  using  a  PC,  as 
illustrated  here,  are  certainly  reasonable. 

The  process  of  calculating  kj  for  particular  dielectric  radii  can  be  automated, 
enabling  kz  to  be  plotted  as  a  function  of  a  continuously  increasing  (decreasing)  radius — 
something  similar  to  a  function  generator.  This  is  quite  useful,  among  other  things,  in 
mode  identification.  As  a  standard,  the  subscripts  on  HE^  and  EH^  refer  to  the  n*  order 
and  the  m*  rank  [18].  In  a  homogeneous  guide,  the  rank  is  determined  by  the  successive 
ordering  of  zeros  for  Jn.  In  the  inhomogeneous  case  this  no  longer  holds.  The  method  of 
identification  used  here,  as  discussed  in  Section  2.5,  is  to  use  the  homogeneous  guide  as  a 
limiting  case.  Thereby,  the  dielectric  radii  within  the  guide  are  gradually  made  larger  and 
the  kz  versus  radius  plot  is  made.  The  starting  point  for  this  plot  is  chosen  as  the  desired 
TE  or  TM  mode  in  the  homogeneous  guide,  usually  air  filled.  This  method  of  tracing  out 
the  1^  has  a  number  of  advantages  over  successive  applications  of  MDCW  including  speed 
and  accuracy.  A  speed  increase  is  quite  obvious  since  all  of  the  parameters  do  not  need  to 
be  entered  for  each  increment  in  radius.  Less  obvious  is  the  greater  speed  realized  by  using 
the  previous  last  two  calculations  for  kj  along  with  a  new  calculated  guess.  This  generally 
halves  the  number  of  iterations  needed  for  the  calculation  of  the  new  root.  Accuracy  is  also 
increased  since  it  is  much  easier  to  stay  with  the  same  mode  when  tracing  than  with 
MDCW.  One  can  accidentally  jump  to  other  modes  with  the  same  order  but  different  rank 
if  the  initial  guesses  are  too  far  from  the  desired  mode. 

3.3  Results  for  Two  Dielectrics 

At  this  point  it  may  be  instructive  to  consider  a  few  examples  and  compare  the 
results  with  published  data.  Two  modes  of  particular  interest  are  the  TEn  and  TMn  in 
addition  to  their  hybrid  counterparts  the  HE,,  and  EH, ,  modes.  In  the  homogeneous  case, 


the  TEn  mode  is  fundamental  and  the  HEU  remains  so  for  the  inhomogeneous  case, 
hence,  its  importance.  The  importance  of  the  EHn  mode  will  become  more  evident  later  in 
this  chapter  and  in  Chapter  5  when  the  experimental  measurements  on  the  finite  guide  are 
discussed.  Consider  a  0.24  m  diameter  pec  waveguide  with  a  dielectric  filling  having 
£rj=l0  and  £r2=l  F/m  as  in  Figure  3.2.  For  the  TEn  mode  in  an  air-filled  homogeneous 
guide,  k^=14.3  rad/m  at  1  GHz.  There  are  two  ways  in  which  a  versus  relative  radius 
plot  can  be  constructed  as  the  insert  dielectric  rod  becomes  larger  (or  smaller)  in  diameter. 
One  way  is  with  successive  applications  of  MDCW  while  gradually  increasing  (decreasing) 
the  radii.  The  best  way  is  through  automation  of  this  process  and  allowing  the  computer 
program  to  perform  these  calculations.  (As  has  already  been  mentioned,  due  to  space 
limitations,  this  program  is  not  listed  in  this  thesis.) 

Starting  the  tracing  process  for  an  exceedingly  small  radius,  the  plot  in  Figure  3.3  is 
constructed.  The  HEU  mode  is  christened  in  the  limit,  as  a/b— >  0  or  1  the  TEn  mode  is 
obtained  for  the  homogeneous  guide.  One  way  this  graph  can  be  interpreted  is  as  a 
transition  diagram  from  one  homogeneous  case  to  another.  The  is  a  function  of  the  inset 
radius  and  is  known  analytically  only  for  the  homogeneous  cases,  i.e.,  at  the  endpoints. 
This  graph  indicates  how  the  axial  wavenumber  varies  between  the  two  homogeneous 
situations.  This  plot  compares  well  with  one  given  by  Harrington  [21]  and  is,  in  fact,  an 
exact  match.  The  transverse  electric  fields  for  the  four  cases  labeled  A,  B,  C  and  D  are 
shown  in  Figure  3.4.  In  each  case,  Ep  is  discontinuous  at  the  boundary  by  the 
accumulation  of  bound  charge  at  this  interface.  The  field  is  particularly  interesting  in 

that  the  field  maximum  begins  outside  the  dielectric  rod  and  gradually  moves  inside  as  the 
rod  radius  is  increased.  Related  to  this  is  the  transformation  of  kp2  from  a  purely  real  to  a 
purely  imaginary  quantity  at  kx/k0=l.  The  wavenumbers  are  listed  in  Table  3.1  for  these 
four  cases.  The  region  much  below  this  point  is  characterized  by  electric  fields  primarily 
concentrated  in  the  outer  dielectric,  whereas  the  region  much  above  kz/ko=l  has  its  electric 
fields  concentrated  in  the  inner  dielectric.  The  magnetic  fields  retain  most  of  their  shape  as 
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a/b  is  increased  with  the  maximum  fields  located  in  the  center  dielectric.  The  transverse 
electric  fields  also  display  an  increasingly  exponential  behavior  in  region  2  as  a/b  increases. 
This  is  a  further  manifestation  of  the  imaginary  behavior  of  kp2  which  gives  rise  to 
"trapped"  waves  in  the  inner  dielectric.  Imaginary  arguments  for  the  Bessel  functions  in 
(2.14)  give  linear  combinations  of  modified  Bessel  functions  which  vary  smoothly  as 
opposed  to  the  oscillatory  behavior  of  the  normal  Bessel  functions.  Behavior  of  this  son  is 
reminiscent  of  surface  waves  and  demonstrates  itself  here  for  kp2  sufficiently  large  and 
imaginary.  If  there  were  a  coating  instead  of  an  insert  (meaning  £ri=l  and  £r2>l),  then  the 
situation  would  be  reversed.  Discussions  of  this  arrangement  are  given  by  Chou  and  Lee 
[7]  and  Lee  era/.  [28]. 

Another  interesting  example  is  an  air-filled  waveguide  of  1  in  diameter  having  a 

dielectric  insert  of  eri=37.6  F/m.  For  the  TEn  mode  in  an  air-filled  guide,  k^==-j  118.3 

rad/m  at  4  GHz  serves  as  a  starting  point  for  the  kz  versus  relative  radius  plot,  with  the 

results  of  this  given  in  Figure  3.5.  This  graph  has  a  peculiar  shape  in  that  the  HEn  and 

EHn  modes  tend  to  form  a  continuous  curv**  This  is  in  stark  contrast  to  these  similar 

plots  for  the  homogeneous  case.  Also,  a  backward-wave  region  (where  the  phase  and 

group  velocities  have  differing  signs)  has  formed  in  the  kz  versus  a/b  plot  where  the  slope 

aco 

is  negative.  This  backward-wave  assignment  can  be  arrived  at  here  by  rewriting  a/b  as 
with  c=speed  of  light  and  a  is  a  constant.  The  x  axis  can  now  be  thought  of  as  an 
increment  in  0)  instead  of  a/b.  Although  this  graph  was  generated  as  a  (3  versus  relative 
radius  plot,  the  x  axis  can  be  regarded  as  proportional  to  frequency  from  which  backward- 
wave  assignments  can  be  made  by  observing  the  slope  so  long  as  the  constitutive 
parameters  are  not  frequency  dependent. 

Another  interesting  feature  of  this  graph  is  the  region  from  a/b=0.302  to  0.358 
where  no  purely  real  or  purely  imaginary  solutions  for  kz  exist.  Instead,  complex  solutions 
are  found  in  the  form  of  (3.2).  These  solutions  are  the  complex  or  "leaky”  modes 
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discussed  in  Chapter  1  and  have  been  investigated  in  some  detail  in  the  literature  [14],  [29], 
[30].  For  conservation  of  energy  reasons,  these  modes  exist  in  pairs  having  negative 
complex  conjugate  propagation  coefficients  k^  and  -k^*.  In  this  trace,  the  backward-wave 
and  complex  mode  regions  exist  simultaneously.  Clarricoats  and  Taylor  [14]  show  that 
this  is  always  the  c?se  for  the  two-dielectric  inhomogeneity  but  the  converse  does  not 
apply.  That  is,  complex  modes  can  exist  without  propagating  backward-wave  regions. 
Later  in  Section  4.2  this  same  occurrence  will  be  shown  to  exist  in  the  three-dielectric 
inhomogeneity.  Since  there  is  no  fundamental  difference  between  three  dielectrics  and  four 
or  more,  it  can  be  safely  stated  that  for  any  number  of  dielectric  layers,  the  existence  of 
propagating  backward  waves  in  a  k^  trace  is  always  accompanied  by  a  complex  mode 
region,  but  not  necessarily  vice  versa. 

Once  this  plot  has  been  constructed,  the  transverse  fields  can  be  computed  using  the 
kz  value  for  the  desired  relative  radius.  Choosing  a/b=0.788,  then  from  Figure  3.5, 
kz=461.7  rad/m  for  the  HEn  mode  and  kj,=265.1  rad/m  for  the  EHn  mode.  The  fields  for 
these  two  cases  as  computed  by  MDCW  are  shown  in  Figures  3.6  and  3.7.  They  are  all 
normalized  so  that  the  maximum  field  component,  either  E  or  H,  has  a  magnitude  of  one. 
Figures  3.6  (a)  and  (b)  show  the  transverse  electric  and  magnetic  fields,  respectively,  for 
the  HEn  mode.  All  of  the  fields  are  continuous  except  for  the  Ep  component  which  is 
discontinuous  by  the  accumulation  of  bound  charges.  Similar  plots  for  the  EH, ,  mode  are 
shown  in  Figure  3.7.  These  graphs  compare  exactly  to  those  of  Zaki  [4]  except  for  a 
different  normalization  of  the  H  fields.  Zaki  also  subscribe?  to  some  nonstandard 
nomenclature  for  the  modal  designations  [3].  Graphs  of  these  mode  patterns  are  shown  in 
Figure  3.8.  Figure  3.6  agrees  with  Figure  3.8(a)  and  Figure  3.7  with  3.8(b).  The 
excellent  correspondence  shown  here  serves  as  a  good  source  of  verification  for  the 
accuracy  and  validity  of  the  results  from  the  programs  for  two  dielectrics.  For  more  than 
two,  verification  will  be  provided  by  experimentation.  It  is  not  surprising  that  such  good 
agreement  has  been  reached.  The  analysis  used  here  and  that  of  Zaki  are  founded  on  the 
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same  approach,  namely,  expression  of  the  eigenmodes  as  a  product  of  the  expansion 
functions  as  in  Chapter  2.  Figure  3.8  also  contains  the  field  patterns  for  the  transverse  E 
and  H  fields  across  the  waveguide  cross  section. 

Up  to  this  point,  fairly  typical  topics  have  been  discussed  concerning  the 
multidielectric  waveguide.  Now,  some  additional  discussions  will  be  given  for  the  two- 
dielectric  case.  These  examples,  which  include  k^  versus  a),  versus  Er,,  and  3-D  power 
flow  plots,  will  lend  further  insight  into  the  rod-insert  guide  and  also  bridge  the  way  for 
three  and  more  dielectrics. 

3.4  Other  Axial  Wavenumber  Plots 

The  first  of  these  plots  to  be  demonstrated  is  the  k^-co  plot  for  the  frequencies  1  to 
10  GHz.  Construction  of  this  begins  by  determining  the  starting  point  for  the  frequency 
and  mode  of  interest.  This  is  done  by  tracing  out  the  k^  versus  relative  radius  plot  at  I 
GHz.  The  results  of  this  are  shown  in  Figure  3.9  for  b=0.5  in,  £^=37.6  and  Er2=l  F/m. 
The  two  modes  shown  arc  HEn  and  EH^.  Again,  a  region  of  complex  mode  solutions 
exists  as  for  the  4  GHz  case  except  that  here  it  connects  two  evanescent  hybrid  mode 
regions.  It  is  interesting  to  observe  that  in  this  complex-mode  region  the  distinction 
between  HEn  and  EHn  disappears  and  both  modes  seem  to  exist  simultaneously.  It  is  in 
this  complex-mode  region  where  the  starting  point  for  the  kz-0)  plot  is  obtained  at 
a/b=0.788.  Since  HEn  ^  EHn  appear  to  exist  simultaneously  at  this  point,  it  suffices  as 
a  starting  point  for  both  k^-to  plots.  The  resulting  kz-co  curves  are  shown  in  Figure  3.10. 
At  a  frequency  slightly  greater  than  2  GHz,  the  modes  revert  from  complex  to  hybrid  with 
the  EHn  experiencing  a  small  backward-wave  region  from  approximately  2  to  2.4  GHz. 

Also  in  Figure  3.10  are  the  J3-co  plots  above  cutoff  for  the  TEn  and  TMn  modes 
in  the  homogeneously  filled  guide  of  1  in  diameter  with  er=37.6  F/m.  At  lower 
frequencies,  one  can  see  a  large  discrepancy  between  the  transverse  and  hybrid  modes  as 
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expected.  At  higher  frequencies,  however,  it  is  quite  surprising  that  the  plots  all  have 
basically  the  same  shape  and  trends.  This  indicates  that  the  inhomogeneous  guide  could  be 
modelled  as  homogeneous  at  frequencies  much  beyond  cutoff,  at  least  for  these  two  cases, 
by  selecting  an  effective  Ej. 

In  addition  to  the  variation  of  kj  as  a  function  of  co,  is  the  important  case  in  which 
kz  varies  as  a  function  of  the  dielectric  Etv  This  is  useful  from  a  number  of  standpoints. 
One,  as  different  loading  materials  are  used,  the  cutoff  frequency  of  the  composite  structure 
can  be  observed.  Two,  the  variation  in  attenuation  afforded  by  each  mode  can  be 
scrutinized  as  the  imaginary  part  of  £r;  (representing  conduction  losses)  is  varied.  An 
example  of  the  former  is  shown  in  Figure  3. 1 1  for  the  two-dielectric  problem  having 
a/b=0.788,  b=0.5  in  and  £r2=l  F/m  at  4  GHz.  The  inner  dielectric's  real  permittivity  is 
varied  from  1  to  38  F/m  for  both  the  HEn  and  EHU  modes.  For  £ri=l  F/m,  both  modes 
revert  to  the  homogeneous  case  which  serves  as  a  convenient  starting  point  for  mode 
classification  reasons.  As  the  dielectric  constant  is  increased,  kj  increases  monotonically  in 
these  two  situations.  Cutoff  occurs  at  Ei-^5.3  and  Eri=13.9  F/m  for  the  HEn  and  EHn 
modes,  respectively. 

3.5  Power  Flow 

Power  flow  is  an  important  topic  in  wavep  aide  considerations  and  will  be  discussed 
here  for  the  multidielectric  case.  In  the  most  general  sense,  the  total  power  across  some 
transverse  plane  along  the  guide  is  given  by 


P  =  J  E  x  H*  •  ds 


(3.6) 


where  the  *  refers  to  complex  conjugation.  In  general,  all  six  components  of  both  E  and  H 
will  exist,  but  since  only  the  z  directed  power  flow  is  of  importance  here 
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and 


S  =  E  x  H* 


v  - 


W 


i=l....,q.  (3.7) 


Considering  Equations  (2.5)  and  (2.6)  for  the  p  components  gives 


kjlcj  „ 

hm(n0)  +  —  vfhm(n<(>) 


(3.8) 


kJc,  . 

■JjJTVi  he(n<J>)  +  -^he(n<D)^ 


(3.9) 


where  the  primes  denote  differentiation  in  accordance  with  (2.19).  Substituting  these  two 
equations  plus  (2.17)  and  (2.18)  into  (3.7)  gives  in  part 
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nkp.*  I  I  2 

o)2ei^*p 

Equation  (3.7)  can  be  used  to  calculate  the  power  flow  in  the  z  direction  by  Equations 
(3.10)  and  (3.1 1).  In  particular,  the  real  part  of  Szj  is  of  primary  importance  in  this  study 
since  it  represents  the  time-averaged  power  flow  in  the  axial  direction. 

An  important  special  case  is  the  azimuthally  symmetric  modes  for  which  n=0. 
Equation  (3.7)  then  reduces  to 


S 
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|kPi[Ci,V(kPip)  +  C£!Nn'(kPip)]|2  . 


(3.12) 


If  the  guide  is  homogeneous,  only  one  coefficient  will  be  left,  namely,  Cn.  To  ensure 
nonzero  power  flow,  this  coefficient  is  chosen  equal  to  one  when  solving  for  the 
coefficients  [C].  This  is  important  since  it  was  found  that  in  the  limit  as  the  guide  becomes 
homogeneous,  choosing  other  coefficients  equal  to  one  may  force  the  Cn  coefficient  to 
zero  in  the  case  where  n=0.  Since  nonzero  Szj  is  to  be  expected,  a  boundary  condition  of 
sorts  is  implied  by  this  seemingly  nonphysical  result. 

An  example  of  the  power  flow  through  a  cross  section  of  a  waveguide  is  shown  in 
Figures  3.12  and  3.13.  These  figures  contain  graphs  of  the  three-dimensional  (3-D)  power 
flow  along  with  a  slice  through  the  middle  of  the  3-D  plot  in  the  0=0  plane.  Figure  3.12 
corresponds  to  the  point  labeled  A  in  Figure  3.3  and  Figure  3.13  corresponds  to  point  D. 
These  two  graphs  show  a  tremendous  redistribution  of  the  power  flow  from  largely  in  the 
outer  dielectric  to  mostly  in  the  inner  dielectric.  For  a/b=0.1,  99.5%  of  the  power  is  in  the 
outer  dielectric  while  for  a/b=0.4,  95%  of  the  power  concentrates  in  the  inner.  Although 
the  radius  of  the  inner  dielectric  increases  fourfold,  at  a/b=0.4,  the  inner  dielectric  still  only 
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comprises  16%  of  the  cross-sectional  area.  A  question  arises  as  to  why  the  power  shifts  in 
concentration.  This  redistribution  of  power  is  marked  by  the  magnitude  and  phase  of  kp2- 
As  kp2  becomes  increasingly  imaginary,  the  power  tends  to  concentrate  in  the  inner 
dielectric  as  noted  by  these  two  figures  and  the  values  in  Table  3.1.  Although  the  kp2=0 
point  at  a/b  =0.244  does  not  mark  a  50/50  power  distribution  between  inner  and  outer 
dielectrics  (this  actually  occurs  at  a/b=0.267),  a  large  imaginary  radial  wavenumber  does 
signal  a  small  real  axial  power  flow  in  that  dielectric. 

Another  example  is  shown  in  Figure  3.14.  This  plot  is  for  the  point  labeled  C  in 
Figure  3.5  having  ^*165  rad/m  in  the  HEn  mode.  Here  a/b=0.36  and  b=0.5  in.  One 
interesting  feature  of  this  power  plot  is  the  existence  of  regions,  in  the  vicinity  of  the  <p=0 
plane,  where  power  is  actually  flowing  backwards  towards  the  source.  In  this  particular 
example,  the  total  net  power  flow  is  positive  with  69%  of  the  power  concentrated  in  the 
inner  dielectric  and  31%  in  the  outer  dielectric  and  negative.  However,  the  inner  dielectric 
covers  only  13%  of  the  cross-sectional  area  of  the  waveguide.  This  is  another  good 
example  of  how  the  axial  power  flow  concentrates  in  the  dielectric  with  a  higher  real  radial 
wavenumber.  This  is  a  rule,  in  general,  but  the  actual  percentages  are  also  a  function  of  the 
mode  type  and  the  relative  dielectric  radii. 

The  negative  power  flow  phenomenon  seen  above  appears  to  be  a  remnant,  so  to 
speak,  from  the  backward-wave  region.  That  is,  from  taking  a  number  of  cases  along  the 
P  versus  relative  radius  plot,  it  was  determined  that  the  power  flow  in  the  center  and  outer 
dielectrics  is  oppositely  directed  for  all  values  of  a/b  in  these  two  modes.  The  results  of 
this  are  shown  in  Table  3.2.  Although  the  last  entry  has  a  positive  value  in  the  outer 
region,  there  still  is  some  negative  flow;  however,  it  is  less  than  the  positive  flow  in  that 
annulus.  The  values  listed  in  Table  3.2  are  ordered  as  one  would  trace  out  an  ascending 
path  along  the  outside  edge  of  the  plot  in  Figure  3.5  beginning  at  point  A.  Examining  these 
percentages,  a  picture  begins  to  emerge  of  the  power  flow  slowly  diminishing  in  the 
backward  direction  within  the  inner  dielectric  and  a  greater  percentage  going  forward. 
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Further  insight  into  this  is  gained  through  consideration  of  Figure  3. 15  which  is  the  power 
flow  plot  for  point  A  in  Figure  3.5  and  the  First  entry  in  Table  3.2.  Here  the  positive  power 
flow  is  in  the  outer  dielectric  and  the  negative  within  the  inner.  Owing  to  the  backward- 
wave  nature  of  the  wave  here,  the  phase  fronts  are  toward  the  source  and  the  energy 
propagation  is  away.  In  the  backward-wave  region,  the  power  flow  is  positive  in  the  outer 
dielectric,  negative  in  the  inner.  The  converse  of  this  applies  to  the  forward-wave  region  of 
the  EHn  mode  at  a/b=0.5815  depicted  in  Figure  3.16. 

This  concept  of  a  "backward  wave"  seems  to  be  somewhat  muddled  for  the  multiple 
dielectric  cases  considered  here  since  power  is  actually  traveling  in  both  directions  in  both 
the  forward-  and  backward-wave  regions!  One  needs  to  keep  in  mind  that  in  the  backward- 
wave  region,  the  phase  and  group  velocities  have  opposite  signs.  To  choose  the  correct 
sign  for  kj,  as  in  Table  3.2,  an  introduction  of  a  small  amount  of  loss  to  the  inner  dielectric 
will  make  the  right  choice  apparent.  Letting  eri=37.6-j0.01  F/m  for  the  case  with 
a/b=0.3646  in  Figure  3.15,  two  roots  to  the  characteristic  equation  were  found  using 
MDCW.  These  are  kz=68.22+j0.179  and  k^-68. 22-jO.  179  rad/m.  The  first  represents  an 
augmented,  forward  traveling  wave  while  the  second  an  augmented  backward  one.  The 
second  is  the  correct  one  since  the  augmentation  is  bounded  by  the  source  and  hence  a 
physical  result.  In  this  backward-wave  region,  the  correct  result  for  is  with  both  p  and 
a  negative.  This  is  the  reason  for  plotting  abs(p)  in  Figures  3.5  and  3.10,  for  example. 

It  is  also  interesting  to  inquire  as  to  which  mode,  either  HEn  or  EHn,  to  name  the 
backward-wave  region.  By  Snitzer's  scheme  for  the  sign  of  Cn/Cl3,  the  mode  is  correctly 
called  EHn  since  the  sign  is  negative.  However,  if  P  were  mistakenly  chosen  as  positive, 
the  HEU  mode  would  result  since  the  sign  of  the  ratio  then  becomes  positive.  In  the  past, 
the  backward-wave  region  would  indeed  be  labelled  HEn.  In  the  next  chapter,  it  will  be 
shown  in  considerable  detail  that  this  region  is  actually  part  of  the  EHj  j  mode  trace. 


Table  3.1. 


Wavenumbers  of  the  four  cases  labelled  in  Figure  3.3  for  the 

HEn  mode  in  the  rod- insert  guide  having  b=0.4X,  f=l  GHz, 

e.  =10  and  e.  =1. 

“  *2 


Case 

kj  (rad/m) 

a 

b 

kpi  (rad/m) 

kp2  (rad/m) 

A 

14.85 

0.1 

64.59 

14.79 

B 

17.39 

0.2 

63.95 

11.70 

C 

34.37 

0.3 

56.66 

-j27.25 

D 

49.36 

0.4 

44.22 

-j44.69 

Table  3.2. 

Axial  wavenumber  and  percent  total  axial  real  power  flow  for  the  five 
cases  labelled  in  Figure  3.5  for  both  the  HEU  and  EHU  modes  in  the 
rod-insert  guide  having  b=0.5  in,  f=4  GHz,  £^=37.6  and  6^=1. 


Case 

kz  (rad/m) 

a 

b 

%  of  total  real  power  flow 

inner 

outer 

A 

-68 

0.3646 

-22.77 

77.23 

B 

-89 

0.3600 

-31.61 

68.39 

C 

165 

0.3600 

68.76 

-31.24 

D 

262 

0.3900 

94.42 

-5.58 

E 

462 

0.7880 

99.78 

0.22 

Figure  3. 1 .  Example  of  the  quadratic  fit  procedure  to  locate 
the  zero  of  a  function  by  Muller's  method 


Figure  3.2.  Rod-insert  waveguide  arrangement. 


Normalized  1^  versus  relative  radius  plot  for  the  HEn  mode 
rod-insert  guide  having  b=0.4  X,  f=l  GHz,  e,  =10  and  ^  =1 
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Figure  3.4.  (b)  Electric  fields  corresponding  to  the  k£  of  point  B  in  Figure  3.3. 


Figure  3.4.  (c)  Electric  Helds  corresponding  to  the  k*  of  point  C  ir  Figure  3.3. 


versus  relative  rod  radius  for  the  HEn  and  EH,,  modes  in 
ert  waveguide  having  b=0.5  in,  f=4  GHz,  e,  =37.6  and  ^  = 


versus  relative  rod  radius  for  the  HEn  and  EHn  modes, 
e  complex-mode  region,  in  the  rod-insert  waveguide  having 


curves  for  both  the  rod-insert  guide  having  a/b=0.788 
and  e,  =1 ,  and  the  homogeneous  guide  of  radius  0.5  in 


varafuld*  llunatcr  (m) 


(a) 


(b) 


Figure  3. 1 5.  Axial  real  power  flow  of  the  EH[  [  mode  (a)  in  the  0=0  plane  and.  (b)  3-D, 
for  point  A  in  the  kj  trace  of  Figure  3.5  having  a/b=0.3646,  b=0.5  in, 
f=4GHz,  e,  =37.6  and  £=1. 

M  “2 


wi>«|«Ui  tltaaUr  (m) 


(a) 


(b) 


Figure  3. 16  .  Axial  real  power  flow  of  the  EHU  mode  (a)  in  the  <}>=0  plane  and,  (b)  3-D 
for  kz=70  in  the  kj  trace  of  Figure  3.5  having  a/b=0.5815,  b=0.5  in, 
f=4  GHz,  £,=37.6  and  ^=1. 
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CHAPTER  4 

INHOMOGENEOUS  THREE-LAYERED  DIELECTRICALLY- 
LOADED  WAVEGUIDES 

The  theoretical  analysis  used  to  formulate  and  construct  the  computer  programs  for 
the  circular  waveguide  filled  with  layered,  concentric  dielectric  media  is  general  enough  to 
find  solutions  for  any  number  of  dielectrics  having  real  or  complex  constitutive  parameters. 
In  the  preceding  chapter,  the  cases  where  two  dielectrics  filled  the  waveguide  were 
compared  with  the  previously  published  work  as  a  verification  of  the  computer  codes  in 
addition  to  investigating  the  wave  behavior  for  the  rod-insert  arrangement.  In  this  chapter, 
the  waveguide  will  be  assumed  to  be  filled  with  three  concentric,  perhaps  lossy,  dielectrics. 
The  three-dielectric  problem  is  a  reasonable  way  in  which  a  waveguide-below-cutoff 
penetration  into  a  shielded  enclosure  can  be  modelled  when  a  hose  filled  with  some  material 
is  passed  through  the  waveguide,  into  the  shelter.  Since  the  three-dielectric  problem  is  so 
important,  a  considerable  amount  of  emphasis  will  be  placed  on  investigating  the  wave 
behavior  in  such  an  environment  and  obtaining  representative  samples  of  the  key 
characteristics,  such  as  attenuation. 

The  amount  of  attenuation  the  wave  undergoes  in  the  waveguide  is  the  most 
important  characteristic  in  shielding  applications.  Hence,  the  wave  attenuation  will  be 
calculated  as  a  function  of  both  relative  dielectric  radius  and  frequency  for  varying 
permittivities.  Only  losses  from  the  dielectrics  will  be  considered  here  since  the  typical 
waveguides  are  short  (6-12  in)  so  the  metal  waveguide  wall  losses  will  be  relatively  small. 
A  very  important  result,  which  will  be  shown  in  Section  4.3,  is  that  if  the  dielectric  losses 
increase  in  the  guide,  the  wave  attenuation  is  not  necessarily  increased  and  may  actually  be 
decreased!  This  remarkable  finding  is  due  to  a  redistribution  of  power  flow  into  less  lossy 
dielectrics,  that  is  not  possible  in  homogeneous  waveguides. 


Also  in  this  chapter,  a  new  modal  designation  will  be  presented  for  propagating 
modes  in  the  backward-wave  regions.  This  new  designation  will  be  shown  to  be  more 
consistent  with  other  properties  of  the  wave,  such  as  net  power  flow  and  boundedness  of 
the  wave.  A  spinoff  from  this  discussion  is  an  investigation  into  perhaps  a  refined 
definition  for  the  modal  cutoff  frequency  when  backward  waves  are  present  in  the 
dispersion  graph. 

4.1  Lossless  Three-Dielectric  Waveguides 

In  most  applications  of  the  waveguide  below  cutoff  where  a  fluid  is  transferred  into 
a  shielded  enclosure,  the  waveguide  can  be  modelled  as  a  three-dielectric  problem,  that  is,  a 
hose  filled  with  some  material  inside  an  air-filled  guide.  Theoretical  investigations  into  this 
arrangement  proceed  almost  identically  to  those  in  Chapter  3  for  the  two-dielectric  guide. 
The  computer  programs  were  modified  slightly  to  reflect  a  hose  as  the  second  dielectric  by 
keeping  the  radius  difference  between  the  inner  dielectric  and  that  of  the  hose's  outer  radius 
constant  as  the  inner  radius  is  increased  in  a  kj  versus  relative  radius  plot.  An  example  of 
the  results  for  the  kz  versus  a/b  tracing  process  in  the  three-dielectric  problem  is  given  in 
Figure  4.1  for  the  lossless  case.  Here  a  0.1908  in  thick  hose  of  £j=3  F/m  is  increased  in 
diameter  inside  a  2.005  in  radius  waveguide  at  1  GHz  while  the  permittivity  of  the  filling 
inside  the  hose  is  varied  for  the  HEU  and  EHU  modes.  All  th.  ;  curves  vary  smoothly 
until  a/b=*0.9  where  the  outer  edge  of  the  hose  meets  the  waveguide  wall.  For  radii  greater 
than  this,  the  hose  is  gradually  "pushed  out"  of  the  guide  until  a/b=l  where  the  guide 
becomes  homogeneous.  This  a/b=l  point  is  a  convenient  check  for  the  correctness  of  the 
mode  classification  as  discussed  in  Section  3.3.  These  versus  relative  radius  plots,  in 
addition  to  their  interpretation  as  k^-m  plots,  can  also  be  viewed  as  transition  diagrams  of 
the  axial  wavenumber  as  the  guide  changes  from  an  inhomogeneous  guide  of  (q-1) 
dielectrics  to  a  homogeneous  guide  filled  with  a  material  having  parameters  (pn,erj). 
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If  the  permittivity  of  the  hose  is  increased  from  3  to  10  F/m  while  all  other 
parameters  are  left  unchanged,  a  new  versus  relative  radius  plot  shown  in  Figure  4.2 
results.  There  is  little  difference  between  this  plot  and  the  one  in  Figure  4. 1  except  when 
eri  reaches  7  F/m  where  the  HEU  and  EHn  mode  traces  form  continuous  curves.  A 
complex-mode  region  develops  between  these  two  curves  from  a/b=0.5  to  0.6. 
Interestingly,  a  complex-mode  region  has  formed  but  the  propagating  backward-wave 
region  usually  associated  with  the  complex  mode  is  absent.  This  has  been  reported 
elsewhere  [14]  and  this  graph  provides  additional  verification  of  the  phenomenon  when 
more  than  two  dielectrics  fill  the  guide. 

The  EHt !  modes  in  both  figures  contain  regions  of  negative  slope  indicative  of  a 
backward-wave  evanescent  mode.  Two  of  these  regions  exist  bracketed  by  forward- wave 
evanescent  mode  regions.  Although  purely  evanescent  here,  the  introduction  of  losses  will 
impart  propagating-wave  behavior  onto  the  fields  and  the  amount  of  attenuation  in  these 
modes  will  then  become  important.  The  HEn  modes  do  not  exhibit  this  backward-wave 
behavior  in  the  evanescent  region,  but  in  Figure  4.2,  for  a/b  greater  than  approximately 
0.9,  propagating  backward  waves  are  evident.  (An  exception  to  this  is  when  en=2  F/m 
and  a/b>0.9,  when  the  HEU  mode  does  have  a  slight  negative  slope.) 

When  the  permittivity  of  the  innermost  dielectric  is  made  larger  than  approximately 
10  F/m,  the  tracing  process  for  the  versus  relative  radius  plots  becomes  complicated  by 
the  interactions  of  modes  with  the  same  azimuthal  variation  but  differing  radially.  Figures 
4.3  through  4.6  form  a  series  of  graphs  pictorially  depicting  these  interactions  in  the  same 
waveguide  arrangement  as  above,  but  en  ranging  from  10  to  30  F/m.  Both  the  HE  and  EH 
modes  with  unity  azimuthal  variation  are  shown  for  the  radial  variations  1,  2  and  3.  As  en 
increases,  a  few  mode  traces  break  closed  contours  with  other  modes  and  form  smooth 
traces  which  cover  the  entire  a/b  range  as  with  the  HEn  and  EH,;  modes  in  Figures  4.3 
and  4.4.  Conversely,  some  modes  that  were  smooth  across  the  entire  range  form  closed 
contours,  as  with  the  HE12  and  EHJ2  modes.  Where  the  slope  becomes  infinite  in  any  of 
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these  traces,  complex  modes  begin  and  continue  with  increasing  a/b  until  either  a  proper 
evanescent  or  propagating  mode  is  reached.  (For  clarity,  the  complex  modes  are  not 
indicated  in  these  plots  and  only  purely  propagating  or  evanescent  modes  are  shown.) 

A  difficulty  arises  in  Figure  4.5  where  it  is  not  obvious  how  the  EHn  mode  can 
begin  at  a/b=0,  and  then  end  as  a  TM12  mode  at  a/b=l.  Physically  this  is  not  possible, 
since  by  the  definition  of  the  mode,  the  number  of  p  variations  is  a  constant  integer.  A 
similar  situation  occurs  for  this  mode  in  Figure  4.4  for  a  smaller  en,  but  the  transition 
occurs  between  two  infinite  slope  regions  characteristic  of  the  complex  modes.  It  is  the 
HE12/EH12  and  HE12/EHn  coupled  complex  modes  which  provide  the  vehicle  for  smooth 
mode  traces  across  the  entire  range  of  a/b.  Reasoning  similar  to  this  explains  how  the 
mode  transition  from  EHn  at  a/b=0  to  TM12  at  a/b=l  can  take  place  in  Figure  4.5  since  it  is 
these  same  complex  modes  which  are  providing  the  bridge  for  transition  from  one  trace  to 
another.  This  reasoning  is  further  justified  by  Figure  4.7  where  the  complex  modes  are 
now  included  in  the  figure.  The  axial  wavenumbers  for  each  mode  should,  and  can,  be 
traced  by  a  smooth  curve  across  the  entire  range  of  increasing  relative  dielectric  radius. 
However,  as  shown  here,  some  complex  modes  must  be  included  occasionally  to  allow  for 
these  smooth  transitions.  Without  them,  seemingly  nonphysical  traces  can  result  which  are 
difficult  to  explain  as  was  seen  in  Figure  4.5. 

By  increasing  erl  further  to  30  F/m  the  traces  in  Figure  4.6  emerge.  The 
complexity  in  the  interactions  among  the  various  modes  increases  a  great  deal  in  this  graph 
from  the  previous  examples.  The  HE12/EH12  closed  contour  in  the  previous  graph  has 
separated  and  now  the  HE]3  and  EH12  modes  form  closed  contours  in  addition  to  the 
HE!  j/EHn  modes.  Without  graphing  these  traces  in  this  figure  and  those  in  the  past  three, 
it  would  be  difficult  to  properly  identify  some  segments  in  the  traces  of  this  figure,  such  as 
the  segment  S.  However,  by  examining  the  other  traces  and  observing  the  progression 
with  increasing  eri  (which  can  be  thought  of  as  "pegging"  the  left-hand  side  of  the  traces 
and  slowly  increasing  the  right-hand  edge)  this  segment  can  be  identified  as  an  HE,3  mode 
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trace.  Although  the  mode  is  evanescent  in  this  region,  for  lossy  materials,  this  mode  and 
those  nearby  in  the  graph  will  propagate  with  some  large  amount  of  attenuation. 
Knowledge  of  the  behavior  of  these  lower-order  modes  is  important  even  though  they  are 
evanescent  in  the  lossless  case. 

4.2  Complex  and  Backward-Wave  Modes 

The  question  arises  as  to  what  mode  the  backward-wave  region  in  Figure  4.5 
(which  includes  the  point  A)  should  be  named.  Following  the  past  work  of  Clamcoats  and 
Waldron  [9],  [13],  this  region  would  be  designated  the  HEn  mode —  in  particular,  a 
"paranormal  sub-mode"  in  keeping  with  Waldron's  denotation  [31].  This  idea  seems 
reasonable  since  as  the  permittivity  of  the  loading  dielectrics  is  increased,  the  HE!  {  mode's 
(3-0)  curve  gradually  changes  from  that  of  a  homogeneous  guide  to  one  containing  a 
backward- wave  region  for  the  inhomogeneous  case  as  in  Figure  3.10  and  [l3j.  Waldron, 
among  others,  expresses  the  belief  that  the  whole  curve  belongs  to  the  HEn  mode,  but  in 
the  region  below  cutoff,  in  the  backward-wave  area,  kz  is  double- valued  for  this  mode,  so 
that  the  EHU  mode  apparently  "disappears"  for  the  range  of  values  where  the  backward- 
wave  region  exists.  It  will  be  worthwhile  at  this  point  to  investigate  whether  or  not  this  is  a 
correct  assessment. 

The  wavenumbers  and  (Cn.Cn)  coefficients  for  a  number  of  points  labelled  in 
Figure  4.5  of  the  HEn  and  EHn  modes  in  the  propagating  region  are  listed  in  Table  4.1. 
In  order  to  choose  the  correct  kz  root  for  these  points  and  especially  point  A,  a  small 
amount  of  loss  is  added  to  the  dielectrics  (en=20-j0.5,  £i2=3-j0.1  F/m).  From  the  table,  it 
is  apparent  that  the  sign  of  the  real  part  of  C13  for  point  A  agrees  with  points  B,  C  and  D. 
all  of  which  belong  to  the  EH^  mode  trace.  In  keeping  with  Snitzer's  scheme  (cf.  Section 
2.5),  the  backward-wave  region  is  labelled  EHU.  If  on  the  other  hand,  the  root  kz  for 
point  A  had  been  chosen  on  the  other  branch,  which  for  the  backward  wave  :.s  kz=9.3+jl.5 
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rad/m,  the  new  C13  coefficient  would  be  40.4+jl.l.  The  sign  of  the  real  part  of  this 
coefficient  would  then  agree  with  the  HEU  mode  but  the  net  power  flow  would  be 
negative.  Since  the  net  power  flow  is  to  be  away  from  the  source,  the  first  root  is  the 
correct  one  (it  also  has  positive  net  power  flow)  and  the  backward-wave  region  is  rightly 
labelled  an  EHn  mode  and  nowhere  in  the  kz  trace  does  the  EHn  mode  disappear.  This 
idea  was  illustrated  earlier  in  Section  3.5  for  the  two-dielectric  guide  in  reference  to  the  sign 
of  the  net  power  flow.  Depending  on  which  branch  k2  was  chosen  to  be  on,  the  3-D 
transverse  power  flow  plot  appeared  similar  to  either  the  HEn  or  EHn  mode  graphs  from 
other  points  on  the  dispersion  plots.  The  correct  choice  was  made  by  examining  the  type  of 
wave  propagation  present  (either  forward  augmented  or  backward  augmented)  in  addition 
to  the  sign  of  the  net  power  flow. 

An  argument  against  this  type  of  modal  designation  is  that  the  HEn  mode  does  not 
appear  to  have  a  cutoff  frequency  which  is  not  to  be  expected  for  a  propagating  mode.  In 

the  usual  lossless  case,  the  cutoff  frequency  is  defined  by  kz=0.  From  the  separation 

to  5o) 

relation  (2.12),  the  phase  velocity,  vp=-^-,  is  then  infinite  and  the  group  velocity,  vg=0p". 


is  zero  at  cutoff.  Since  the  group  velocity  can  be  a  measure  of  energy  velocity  under  the 
appropriate  conditions,  the  criterion  of  vg=0  may  be  more  fundamental  than  vp=°°  for 
defining  the  cutoff  frequency.  In  Figure  4.7  for  example,  the  HEn  mode  has  vg=0  but 
vp?*0  at  a/b=0.48  where  the  slope  becomes  infinite.  This  point  could  represent  a  cutoff 
frequency  of  the  two  coupled  complex  modes  as  a/b  increases  and  the  cutin  frequency  of 
two  other  modes,  the  forward- wave  HEn  and  the  backward- wave  EHU. 

The  plausibility  of  such  a  notion  as  defining  the  cutoff  frequency  by  vg=0  and  not 
necessarily  with  vp=oo  can  be  investigated  simply  by  considering  the  separation  relation  in 
the  definition  of  vg  as 
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After  some  manipulation  and  noting  that  Re{z} 


z+z* 

2 


4|kJ2 


3kJ  f  3kP  ) 

2a>HiEi(kz+kz*)  -  k,*-^  -  k^-j^J 


(4.2) 


By  defining  the  cutoff  frequency  as 
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(4.3)  is  satisfied,  by  considering  (4.2),  in  the  two  cases  when 

(i) 

|kj2  =  0  =>  p2  =  -  a2 

=>  P  =  a  =  0  for  a,p  real  (4.4) 

and 


(ii) 


— >  00 


(4.5) 


Equation  (4.4)  is  the  usual  definition  for  the  cutoff  frequency,  kz=0.  Equation  (4.5), 
however,  is  a  new  relation  valid  only  in  dispersive  situations,  since  for  a  homogeneous 
guide,  kp  is  a  constant  depending  only  on  the  "electrical  diameter"  of  the  waveguide. 
Therefore,  in  a  homogeneous  guide,  (4.5)  will  never  be  satisfied.  The  kpi  becomes 
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frequency  dependent  in  the  inhomogeneous  guide  because  kz  is  now  the  known  variable 

from  which  kpi  is  calculated  and  not  the  other  way  around  as  in  the  homogeneous  guide. 

Clearly,  kz  is  frequency  dependent  as  any  kz-co  plot  of  the  inhomogeneous  guide  will 

9k* 

attest.  In  an  inhomogeneous  guide,  (4.5)  will  be  satisfied  and  vg=0  whenever  =  «*. 

Although  this  was  known  beforehand.  Equation  (4.2)  lends  some  insight  into  why  an 
alternative  definition  of  the  cutoff  frequency  is  necessary  in  the  inhomogeneous  waveguide. 
In  this  light  then,  the  homogeneous  guide  criterion  for  cutoff,  vg=0  and  vp=°o,  is  a  special 
case  of  the  more  general  situation  where  cutoff  is  defined  by  vg=0  and  vp**»  in  the 
inhomogeneous  guide. 

Additional  confirmation  for  this  method  of  naming  the  ba  1-wave  modes  is 
provided  by  the  graphs  in  Figures  4.8  through  4.1 1.  Here  the  HEn  and  EHU  mode  plots 
are  given  for  the  three-dielectric  inhomogeneity  with  a  small  amount  of  loss  as  e,/  is  varied 

from  10  to  30  and  £^"=0.5  F/m  for  The  physical  parameters  used  to  generate 

these  graphs  are  the  same  as  before  except  for  the  hose,  which  has  a  new  thickness  of  0.2 
in,  that  will  also  be  used  in  the  next  section  for  attenuation  calculations  in  the  three- 
dielectric  lossy  guide.  This  thickness  is  close  enough  to  the  previous  one  of  0.1908  in  for 
an  accurate  comparison.  All  these  plots  illustrate  the  splitting  of  the  complex  modes  with 
the  introduction  of  loss  into  the  dielectrics.  In  Figure  4.8,  the  modes  progress  smoothly 
from  evanescent  to  complex  and  then  either  to  propagating  for  HEn  or  back  to  mostly 
evanescent  for  the  EHn  mode.  No  propagating  backward  wave  is  encountered  by  these 
two  modes.  In  the  complex-mode  region,  the  HE n  mode  is  an  augmented  backward  wave 
while  the  EHu  is  an  attenuated  forward  wave.  The  introduction  of  loss  produces  some 
phase  difference  between  the  two  complex  waves  such  that  nonzero  net  power  flow  is 
possible. 

The  appearance  of  the  backward- wave  region  begins  in  Figure  4.9  for  en=15-j0.5 
F/m.  The  EHn  mode  becomes  backward  augmented  here  in  contrast  to  the  previous  figure 


when  the  HEn  mode  was.  Indeed,  the  backward-wave  region  in  Figure  4.9  has  p 
negative  as  previously  predicted.  This  appears  even  more  reasonable  here  since  a  smooth 
trace  for  p,  from  a  negative  value  in  the  complex-mode  region  for  conservation  of  energy 
reasons  (in  the  limit  as  the  losses  go  to  zero),  to  the  positive  values  for  forward 
propagation,  is  provided  by  negative  values  of  P  in  the  backward-wave  region.  Figures 
4.10  and  4.11  contain  plots  displaying  similar  behavior  for  en=20-j0.5  and  erl=30-j0.5 
F/m,  respectively. 

In  the  Figures  4.8  through  4.11,  the  backward- wave  propagating  region  only 
appears  when  the  EHU  mode  changed  from  evanescent  to  propagating  as  a/b  is  varied  from 
zero  to  one.  As  was  mentioned  in  Section  3.3,  if  anywhere  in  the  k,  trace  there  exists 
propagating  backward  waves,  this  necessarily  implies  that  somewhere  else  in  the  trace 
complex  modes  exist.  Considering  these  traces  then,  it  is  evident  that  the  backward-wave 
region  is  needed  to  preserve  the  continuity  in  the  p  trace  and  assure  that  conservation  of 
energy  is  satisfied  in  addition  to  the  boundedness  of  the  wave.  In  summary  of  this,  under 
certain  conditions,  complex-wave  regions  in  the  trace  can  exist  without  backward-wave 
regions  appearing.  However,  if  backward-wave  regions  do  exist,  somewhere  in  the  trace 
complex  modes  form  since  there  is  no  other  way  for  the  mode  to  end  as  evanescent,  for 
decreasing  a/b,  and  still  preserve  continuity  in  the  k,,  trace. 

One  very  curious  phenomenon  apparent  in  the  plots  having  a  backward-wave 
region  is  the  metamorphosis  of  the  modes,  as  a/b  is  increased,  from  either  evanescent  HEU 
or  EHU  modes,  to  not  propagating  HEU  and  EHn  modes  as  would  be  expected,  but 
rather  to  propagating  EHn  or  HEn  modes,  respectively.  This  behavior  is  not  witnessed  in 
Figure  4.8  where  there  is  no  propagating  backward-wave  region.  As  an  example  of  this, 
starting  with  the  HEn  mode  at  kL=-]ll  rad/m  in  Figure  4.9,  as  a/b  is  increased  from  zero 
to  one,  the  mode  is  first  evanescent  HEu,  then  backward  augmented  in  the  complex-mode 
region.  The  mode  next  becomes  a  propagating  backward  EHj ,  mode  and  finally  a  forward 
propagating  and  attenuated  EHn.  The  gradual  change  from  an  HEn-like  to  an  EHirlike 
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mode  is  clearly  evident  in  this  example.  In  this  line  of  thought  then,  naming  the  backward- 
wave  region  EHn  may  seem  up  to  the  designator  and  not  really  dependent  on  any  physical 
reasoning.  However,  in  addition  to  the  material  already  presented  in  support  of  naming 
this  mode  EH11(  this  trace  clearly  has  propagating  behavior  exclusively  associated  with  the 
EHn  mode,  whereas  the  HEn  mode  is  cut  off  on  this  trace.  The  HEn  mode  displays 
propagating  behavior  only  on  the  other  trace  which  begins  an  evanescent  EHn  in  Figure 
4.9  and  ends  a  propagating  HEU  mode. 

Another  curiosity,  which  will  be  mentioned  in  passing,  is  the  EHn  trace  in  Figure 
4.10.  In  Figures  4.9  and  4.11,  the  EH n  trace  is  "bent”  as  it  passes  from  (3  negative  to 
positive  with  the  associated  bulge  in  a.  In  the  lossless  cases,  as  in  Figures  4.4  and  4.6, 
this  bulge  provides  the  transition  from  a  backward  propagating  EHn  mode  to  a  forward 
propagating  EHn  mode  without  vg  becoming  infinite,  which  would  not  be  physical 
When  this  type  of  behavior  in  the  EHU  mode  is  observed,  as  in  Figure  4.10,  it  has  been 
reported  by  Waldron  [13]  that  a  maximum  backward-wave  region  has  been  attained.  For 
all  other  values  of  eri,  with  i=l,..,q,  the  bulge  will  be  smaller  as  evidenced  by  the  plots  in 
Figures  4.9  and  4.11.  This  phenomenon  has  been  previously  associated  with  the 
HEu/EHn  degeneracy  at  the  point  where  (3=0  as  in  Figure  4.10  [10],  [13].  From  the 
work  that  has  been  presented  here,  the  mode  is  better  named  EHn,  so  there  is  no  mode 
degeneracy  as  both  the  backward-  and  forward- wave  regions  are  EHn.  This  is  certainly 
only  a  modal  nomenclature  and  in  no  way  is  meant  to  invalidate  the  results  of  the  previous 
authors.  The  material  presented  in  this  section  and  in  the  preparatory  ones  is  meant  simply 
to  report  a  new  way  to  name  the  backward-wave  region  that  is  more  consistent  with  other 
properties  of  the  wave  and  modal  designations. 
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4.3  Attenuation  in  Lossy  Three-Dielectric  Waveguides 

The  most  important  characteristic  in  this  study  of  the  three-dielectric  waveguide  is 
the  attenuation  afforded  by  dielectric  losses  above  cutoff.  Once  the  modes  have  been 
properly  identified,  the  losses  for  a  particular  mode  can  be  calculated  as  a  function  of 
relative  radius  or  frequency  with  the  loading  dielectric's  er  as  the  parameter.  The  physical 
dimensions  of  the  waveguide  and  hose,  for  these  theoretical  calculations,  were  chosen  as 
representative  of  a  typical  application  for  the  waveguide-below-cutoff  penetration  in  a 
shielded  enclosure  discussed  in  Chapter  1.  A  4  in  waveguide  of  inner  radius  2.005  in  was 
chosen  and  a  1  in  hose  of  inner  radius  0.4909  in  and  thickness  0.2  in  will  enclose  the 
filling  material  in  the  center. 

The  er  of  the  hose  was  selected  by  consideration  of  the  traces  in  Figure  4. 12.  Here, 
a  is  plotted  against  frequency  for  the  HEn  mode  inside  the  4  in  waveguide  with  the  1  in 

hose  filled  with  a  three-point  dispersive  methanol  having  values  given  in  Section  5.4.  The 
parameter  that  is  varied  among  these  plots  is  the  loss  in  the  hose.  A  value  of  £^'=2.7  F/m 

was  chosen  from  listings  of  values  for  rubber  derivatives  in  von  Hippel  [22].  From  the 
traces  in  Figure  4.12,  a  value  of  £^"=2  F/m  was  selected  since  this  value  placed  the  a -to 

trace  in  with  a  "family"  of  traces  with  higher  z?  's. 

It  appears  from  this  graph  quite  strange  that  more  hose  loss  causes  less  wave 
attenuation.  The  reason  this  occurs  is  as increases,  the  power  flow  is  "expelled"  from 

the  inner  two  dielectrics  and  out  to  the  third,  where  there  is  no  loss  and  the  wave  is 
attenuated  less.  If  the  guide  walls  were  lossy,  this  effect  would  be  less  pronounced; 
however,  in  this  study,  the  guide  is  very  short  so  the  wall  losses  will  be  negligible 
compared  to  the  dielectric  losses.  In  this  section  then,  the  attenuation  afforded  by  dielectric 
losses  will  be  calculated  using  a  value  of  er2=2.7-j2  F/m  for  the  hose  and  all  the 
permittivities  of  the  dielectrics  inside  the  guide  assumed  to  be  independent  of  frequency. 


The  modes  that  will  be  considered  for  these  attenuation  plots  are  the  three  lowest- 
ordered  nonazimuthally  symmetric  modes —  the  HEn,EHu  and  HI^.  Figures  4.13 
through  4.18  contain  traces  of  a  versus  a/b  for  these  three  modes  with  '=  2  and  10  F/m. 

These  two  values  were  chosen  heuristically  and  yield  enough  disparity  to  discern  trends  in 
the  attenuation  as  various  parameters  are  changed.  For  each  e^'  value,  five  values  (1, 

3,  5,  7  and  10  F/m)  were  chosen  as  parameters  for  the  loss.  (Although  these  plots  give  a 

in  Np/m,  conversion  to  dB/m  is  possible  by  multiplying  a  by  8.6859=20  log10(e).) 
Probably  the  most  noticeable  characteristic  of  these  plots  is  that  a  higher  e,.  "  value  does  not 

necessarily  imply  more  loss  as  would  be  the  case  in  a  homogeneous  guide.  These  graphs 
dearly  show  that  the  relative  amount  of  attenuation  between  two  iossy  filling  materials  is 
not  only  a  function  of  the  dielectric  loss  but  also  on  how  large  the  inner  diameter  of  the  pipe 
is  in  relation  to  the  sheath  size.  A  specific  example  of  this  is  in  Figure  4.13  with  a/b=0.6. 
At  point  A,  1^=1 1. 1-j  14. 1  for  eri=2-jl0,  and  at  point  B,  kz=3.7-j24  rad/m  for  en=2-jl 
F/m.  More  attenuation  is  provided  by  the  ?tl=2-jl  F/m  dielectric  than  by  the  other.  This 

relation  reverses  itself  as  the  guide  becomes  homogeneous,  as  it  should,  since  in  this  case  a 
higher  e/'  implies  higher  attenuation  from  dielectric  losses  [20].  Although  only  the  HEU 

mode  was  considered  in  more  detail  here,  all  of  these  modes  in  Figures  4.13  through  4.18 
contain  elements  of  the  same  basic  behavior. 

Another  common  characteristic  among  these  modes  is  the  existence  of  a  region  from 
a/b=0  to  some  point  nearby  where  all  traces  have  approximately  the  same  values  for  a, 
after  which  the  traces  separate.  This  region  is  the  cutoff  region  where  the  modes  are 
propagating  only  slightly  but  are  highly  attenuated.  The  wave  impedance  in  this  region  is 
mostly  reactive  and  the  attenuation  is  primarily  caused  by  the  wavelength  of  the  excitation 
field  being  too  large  to  create  transverse  standing  waves  in  the  waveguide.  Beyond  this 
cutoff  region,  the  attenuation  is  due  to  losses  in  the  loading  dielectrics.  In  a  lossless  guide, 
this  transition  is  abrupt  and  the  cutoff  frequency  is  prominent  at  1^=9.  When  losses  are 
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present,  the  cutoff  becomes  "smeared"  somewhat  and  it  becomes  difficult  to  locate  the 
specific  point  of  cutoff. 

The  other  range  against  which  the  attenuation  will  be  plotted  is  frequency.  Figures 
4.19  through  4.24  contain  plots  of  a  versus  frequency  for  the  same  physical  arrangement 
as  in  the  previous  a  versus  a/b  plots —  specifically,  the  1  in  hose  of  £^-2.7^2  F/tr.  in  the  4 
in  guide.  The  erl  values  used  as  parameters  are  the  same  for  this  graph  series  as  in  the 
previous,  which  are  listed  on  the  plots.  The  frequency  band  is  slightly  different  for  each 
mode  and  was  selected  as  a  good  example  of  attenuation  near  the  cutin  for  each  mode. 

A  common  characteristic  throughout  these  plots,  and  also  very  similar  to  the  a 
versus  a/b  plots,  is  that  the  relative  amount  of  attenuation  between  two  e/'  values  with 
similar  physical  dimensions  is  dependent  not  only  on  the  sizes  of  the  zf  ’s  but  also  on  the 
frequency.  This  is  equivalent  to  the  conclusion  reached  with  the  a  versus  a/b  plots 
provided  er  is  not  a  function  of  frequency,  which  it  is  not  here.  A  similar  statement  can  be 
made  for  the  relation  between  zf  and  the  amount  of  attenuation:  a  larger  e/'  does  not 
automatically  imply  higher  attenuation  for  all  frequencies  in  these  inhomogeneous  guides. 

An  interesting  phenomenon  displayed  in  Figures  4.19,  4.20  and  4.22  is  that  the 
largest  amount  of  attenuation  is  obtained  by  the  smallest  e^"  value.  That  is,  for  £^"=1  F/m, 

much  more  attenuation  is  realized  at  higher  frequencies  than  for  a  larger  £,/'.  The  reason 

for  this  was  discussed  earlier  in  conjunction  with  Figure  4.12,  where  for  the  lo\v\..  ^  " 

value  the  power  flow  is  concentrated  primarily  in  the  inner  dielectric  and  moves  primarily  to 
the  outermost  dielectric  for  larger  values  where  there  is  no  loss.  This  redistribution  of 

power  flow  is  marked  by  the  magnitudes  and  phases  of  the  radial  wavenumbers  in  the 
dielectrics.  As  discussed  earlier  in  Section  3.5,  the  axial  power  flow  concentrates  in  the 
dielectric  with  a  higher  real  radial  wavenumber.  With  losses  present,  the  kpi's  are,  in 
general,  neither  purely  real  nor  purely  imaginary.  Instead,  the  axial  power  flew 
concentrates  in  the  dielectric  having  a  kp,  that  is  large  and  has  a  real  pan  that  is  much  larger 
than  its  imaginary  pan. 
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4.4  Power  Flow  in  Lossy  Three-Dielectric  Waveguides 

It  was  pointed  out  in  the  last  section  that  the  amount  of  loss  present  in  the  loading 
dielectrics  can  have  a  very  large  effect  on  the  power-flow  distribution  in  the 
inhomogeneous  guide.  This  redistribution  can  greatly  change  the  amount  of  attenuation  the 
wave  will  undergo  when  the  dielectrics  within  the  guide  have  largely  differing 
permittivities.  The  effect  that  the  losses  in  the  dielectrics  have  on  the  power  flow  and  wave 
attenuation  can  be  quite  pronounced  depending  on  the  many  factors  mentioned  in  the 
previous  section,  including  the  relative  radius  of  the  inner  dielectrics  in  relation  to  the 
sheath  diameter.  A  good  example  of  this  can  be  illustrated  by  considering  the  four  cases 
marked  in  Figure  4.13.  The  percentages  of  the  total  absolute  real  power  flow  in  each 
dielectric  are  listed  in  Table  4.2.  At  a/b=0.6,  the  power  flow  is  distributed  in  roughly  a 
similar  manner  for  both  cases  A  and  B  having  eri=2-jl0  and  en=2-jl  F/m,  respectively. 
This  is  reaffirmed  by  the  3-D  plots  of  the  power  flow  in  Figures  4.25  and  4.26  for  these 
two  cases.  However,  the  wave  attenuation  in  each  case  is  largely  different  with  a=14.1 
Np/m  in  case  A  and  24  Np/m  in  case  B,  even  though  the  loss  is  ten  times  less  in  the  second 
case.  When  a/b  is  increased  to  0.88,  these  distributions  change  drastically  and  display  a 
couple  of  interesting  attributes.  For  one,  in  case  C  (en=2-jl  F/m,  a=22.5  Np/m)  nearly 
90%  of  the  total  absolute  real  power  flow  is  in  the  innermost  dielectric  while  in  case  D 
(erl=2-jl0  F/m,  a=33.4  Np/m)  there  is  only  25%.  Secondly,  in  case  D,  nearly  half  the 
power  is  flowing  in  the  pipe  when  it  co\  :rs  only  19%  of  the  cross-sectional  area,  with  the 

I 

remainder  roughly  distributed  between  the  two  other  dielectrics.  The  higher  dielectric  loss 
in  region  one  has  caused  the  power  flow  to  concentrate  more  in  the  pipe  (region  two)  and 
also  in  the  air-filled  region  (region  three).  This  tremendous  change  in  the  distribution  of  the 
power  flow  is  illustrated  in  the  3-D  plots  for  cases  C  and  D  in  Figures  4.27  and  4.28, 
respectively. 
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The  attenuation  characteristics  have  reversed  from  those  at  a/b=0.6,  such  that  now, 
for  a/b=0.88,  the  waveguide  arrangement  having  the  larger  dielectric  loss  also  has  the 
larger  wave  attenuation.  That  is,  case  C  has  a=22.5  Np/m  while  case  D,  which  has  more 
dielectric  loss  in  region  one,  has  a=33.4  Np/m.  By  increasing  the  losses  in  region  one, 
less  power  flows  in  this  region,  in  comparison  to  a  similar  guide  of  lesser  dielectric  loss  in 
region  one  (as  cases  A  to  B  and  D  to  C),  but  the  amount  of  attenuation  is  dependent  on  the 
relative  sizes  of  the  dielectric  diameters  to  the  sheath  size. 

The  effects  that  the  relative  sizes  of  the  dielectric  constants  have  on  the  power- flow 

distribution  and  attenuation  can  also  be  quite  pronounced.  The  four  cases  marked  in  Figure 
4.14  correspond  to  the  four  cases  in  Figure  4.13,  except  e^'  has  been  increased  from  2  to 

10  F/m.  Between  cases  B  and  F  at  a/b=0.6,  the  higher  '  value  unexpectedly  gives  rise  to 

a  smaller  attenuation —  from  24  to  16.9  Np/m.  When  the  radius  is  increased  to  a/b=0.88  in 
cases  D  and  H,  the  wave  is  attenuated  almost  equally  in  both  instances  with  a=  33.4  and 
31.5  Np/m,  respectively.  Furthermore,  cases  C  and  G  show  the  greatest  amount  of 

difference  in  the  power-flow  distributions,  but  surprisingly  have  nearly  the  same 
attenuation.  The  reason  for  this  is  that  both  points  have  £^"=1  F/m,  which  is  near  that  of 

£^"=2  F/m  producing  roughly  the  same  amount  of  loss  in  both  dielectrics  and,  in  addition, 

nearly  the  same  percentage  of  power  flows  in  the  inner  two  dielectrics  in  both  cases  C  and 
G.  These  facts  in  combination  produce  the  effect  of  attenuation  which  varies  only 
moderately  when  the  dielectric  constant  in  region  one  is  increased  by  a  factor  of  five. 
These  three  comparisons  demonstrate  that  the  differences  in  the  dielectric  constants  among 
the  regions  in  the  waveguide  can  have  a  large  impact  on  the  power-flow  distribution  and 
wave  attenuation.  More  importantly,  these  comparisons  also  show  that  there  is  no  simple 
relation  between  the  attenuation  and  power-flow  distributions  to  the  relative  dielectric 
diameters  and  dielectric  constants  for  the  inhomogeneous  waveguide.  The  relationships  are 
very  complicated  and  must  be  obtained  for  each  waveguide  arrangement  individually  as 
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was  done  for  these  examples —  by  tracing  out  the  kz  versus  a/b  plots  for  many  cases  while 
varying  the  physical  parameters. 

Furthermore,  a  large  redistribution  in  power  flow  by  itself  does  not  necessarily 
imply  a  significant  change  in  the  attenuation.  For  between  cases  B  and  C  there  is  a  huge 
change  in  the  power-flow  distribution  with  only  a  small  change  in  attenuation,  while  in 
cases  A  and  D  there  is  a  large  redistribution  of  power  flow  and  a  large  change  in 
attenuation.  A  large  variation  in  the  losses  among  the  dielectrics  is  an  additional  factor 
which  must  also  be  considered  when  making  attenuation  evaluations  for  varying  physical 
parameters  such  as  the  dielectric  radii.  Again  it  should  be  emphasized  that  the  power-flow 
distributions  and  the  wave  attenuation  inside  the  inhomogeneous,  layered  waveguide  are 
very  complicated  functions  of  the  dielectric  radii,  constitutive  parameters,  sheath  size  and 
frequency.  No  simple  rules-of-thumb  are  applicable  to  these  waveguides  and  the  correct 
relationships  must  be  obtained  numerically  by  tracing  out  the  axial  wavenumber  and 
computing  the  desired  fielded  characteristic. 


Table  4.1. 

Axial  wavenumbers  and  field  expansion  coefficients  for  the  seven 
cases  labelled  in  Figure  4.5  of  the  HE2 1  and  EHn  modes  in  the 
lossy  three-dielectric  guide  having  e^iO-jO.S,  £^=3^0.1  and  6^=1 


Case 

(rad/m) 

Qi 

C,3 

m 

-9.3-j  1 .5 

1.0 

-40.4-j  1. 1 

B 

4.9-j0.7 

1.0 

-54.6+j5.5 

C 

23.0-jl.0 

1.0 

-58.0+j2.6 

D 

28.8-j  1 .2 

1.0 

-47.5+j2.0 

E 

42.3-j2.9 

1.0 

33.3+j0.5 

F 

62. 1  -j  1 .8 

1.0 

40.0+j0.4 

G 

78.1-jl.4 

1.0 

66.38+jl.l 

Table  4.2. 


Traces  of  kz  versus  a/b  for  the  HEn  and  EHn  modes  in  the  three- 
dielectric  waveguide  at  1  GHz  with  e,  as  parameter,  =10  and  ^  = 


Traces  of  kz  versus  a/b  for  the  HEIm  and  EHlm  modes  in  the  three 
dielectric  waveguide  with  m=l,  2  and  3,  having  f=l  GHz,  b=2.00 


Traces  of  versus  a/b  for  the  HElm  and  EHlm  modes  in  the  three 
dielectric  waveguide  with  m=l,  2  and  3,  having  f=l  GHz,  b=2.00 
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Traces  of  kz  versus  a/b  for  the  HE,m  and  EHlm  modes  in  the  three 
dielectric  waveguide  with  m=l,  2  and  3,  having  f=l  GHz,  b=2.00 


Traces  of  versus  a/b  for  ihe  I IE)m  and  EH|m  modes  in  the  three 
dielectric  waveguide,  including  complex  modes,  with  m=l,  2  and 
having  f=l  GHz,  b=2.(X)5  in,  £.  -  20,  e.  =  3  and  er  =1 . 


Traces  of  kz  versus  a/b  for  the  HEjj  and  EHU  modes  in  the  lossy  three 
dielectric  waveguide  having  f=l  GHz,  b=2.005  in,  t-=0.2  in,  Er  =  10-j0 
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Traces  of  kz  versus  a/b  for  the  HE, ,  and  EH, ,  modes  in  the  lossy  three 
dielectric  waveguide  having  f=l  GHz,  b=2.005  in,  t=0.2  in,  e,  =  20-j0 
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Figure  4.18.  Traces  of  a  versus  a/b  for  the  IIE^,  mode  in  the  three-dielectric  guide 
having  f=2.5  GHz,  b=2.005  in,  t=0.2  in,  £,'=10,  £^=2.7^2  and  £^=1 

with  £,"  as  parameter. 
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Figure  4.20.  Traces  of  a  versus  frequency  for  the  HE,,  mode  in  the  three-dielectric 
guide  having  a-0.4909  in,  b=2.005  in,  t=0.2  in,  £^'=10,  ^=2.7^2  and 

£^=1,  with  e,  "  as  parameter. 


versus  frequency  for  the  EH|j  mode  in  the  three-dielectric 
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Figure  4.27.  Axial  real  power  flow  of  the  HEjj  mode  (a)  in  the  0=0  plane  and,  (b)  3-D, 
for  point  C  in  the  kz  trace  of  Figure  4.13  having  a/b=0.88,  b=2.005  in, 
f=I  GHz,  Cr  -2-jl,  eT2=2.7-j2  and  6^=1. 
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Figure  4.28.  Axial  real  power  flow  of  the  HEn  mode  (a)  in  the  t|>=u  plane  and,  t’o)  3-D, 
for  point  D  in  the  trace  of  Figure  4. 13  having  a/b=0.88,  b=2.005  in. 


f=l  GHz,  eri=2-jl0,  Zt=2J-]2  and  £^=1. 
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CHAPTER  5 

EXPERIMENTAL  INVESTIGATIONS 

The  subject  matter  up  to  this  point  has  been  primarily  directed  toward  the  theoretical 
and  numerical  aspects  of  the  inhomogeneous  circular  waveguide  and  developing  an 
understanding  of  the  principles  involved.  The  emphasis  now  will  be  shifted  toward  the 
more  physical  concerns  of  this  problem  and  the  experimental  investigations.  The  purpose 
of  this  laboratory  work  is  to  see  how  well  the  theoretical  models  and  methods  used  in  this 
development  match  the  physical  phenomenon.  This  also  provides  a  verification  for  the 
accuracy  of  the  computer  programs  developed  and  written  for  this  project  in  addition  to 
their  utility. 

Two  very  different  experimental  techniques  were  devised  and  used  in  these 
measurements.  One  is  a  resonant  cavity  technique  which  proved  very  useful  for 
investigations  with  near  lossless  materials.  The  second  technique  is  a  finite  open-ended 
waveguide  method  which  is  much  harder  to  model  theoretically  but  is  useful  for  lossy  and 
very  lossy  materials.  A  major  problem  encountered  in  analyzing  these  experiments  was  not 
having  accurate  values  for  the  permittivities  of  the  loading  materials.  This  made  analytical 
verification  difficult  but  not  impossible  as  additional  boundary  conditions  which  appeal  to 
logical  argument  made  predicted  permittivities  for  a  select  number  of  materials  reasonable. 

The  finite  open-ended  waveguide  measurements  have  an  additional  appeal  in  that 
they  are  also  used  in  the  shielding  measurements  discussed  in  Chapter  1 .  In  Section  5.4, 
the  results  of  the  attenuation  measurements  for  this  arrangement  will  be  seen  well  modelled 
by  infinite  waveguides  of  identical  cross  section.  This  is  encouraging  for  further 
justification  of  the  use  of  these  programs.  However,  the  use  of  these  codes  is  limited  to 
explaining  measured  results.  Their  use  as  predictive  tools  is  highly  questionable. 
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5.1  Resonant  Cavity  Method  Justification 

One  technique  for  measuring  the  properties  of  the  dielectncally-loaded  infinite 
waveguide  is  with  a  resonant  cavity  [11].  At  first  glance,  the  finite,  bounded  geometry  of 
the  cavity  may  seem  far  from  that  of  the  infinite  waveguide.  However,  the  inclusion  of  end 
plates  on  the  waveguide  does  not  affect  the  separability  of  the  geometry  inside  the  cavity. 
That  is,  the  interior  of  the  cavity  remains  a  separable  geometry  and  the  interior  fields  can  be 
found  by  using  separation  of  variables  as  in  Chapter  2.  Since  one  would  now  expect  the 
waves  to  be  traveling  back  and  forth  in  the  cavity  and  not  propagating  in  only  one  direction, 
the  choice  of  e-jklZ  for  the  z  variation  must  be  replaced.  The  choice  of 

hOc^z)  =  sin(kjz)  or  cos(kzz)  (5.1) 

allows  for  all  the  transverse  electric  fields  to  meet  the  zero  field  boundary  condition  at  both 
pec  end  plates.  Therefore,  for  zero  tangential  fields  (Ep  and  E^  at  z-  0  and  L,  the  length  of 

the  cavity,  from  Equations  (2.5)  and  (2.6) 

sin(kjL)  =  0  => 

^  =  ^7  p  =  04:14:2,...  .  (5.2) 

In  order  to  satisfy  this  boundedness  criterion  in  the  z  direction,  another  free  parameter 
(either  to,  p  or  e)  must  be  fixed.  Using  (2.12),  assuming  p  and  e  are  known,  the  resonant 
frequencies  can  be  fixed  for  each  p,  thus  enforcing  (5.2).  This  is  in  contrast  to  the 
waveguide  problem  where  the  frequency  is  given  and  the  kz  determined.  This  implies  that 


solutions  to  the  resonant  cavity  problem  must  be  found  in  a  different  manner  than  for  a 
waveguide  of  similar  cross  section. 

With  the  geometry  separable,  the  dispersion  relation  given  in  Equation  (2.12) 
remains  valid  for  the  interior  of  the  cavity.  Plots  of  the  dispersion  relations  for  both  the 
homogeneous  and  inhomogeneous  waveguides  considered  earlier  are  contained  in  Figure 
3.10.  The  methods  for  determining  these  plots  in  those  two  cases  are  dramatically 
different,  but,  nonetheless,  they  are  both  graphical  depictions  of  the  dispersion  relations.  If 
metal  end  caps  are  placed  perpendicular  to  the  guide  axis  at  any  two  points  along  z,  a 
resonant  cavity  is  formed  for  which  the  fields  still  obey  the  same  dispersion  relation  as  the 
guide  Fields.  The  additional  boundary  condition  in  (5.2)  forces  a  discretization  of  these 
dispersion  graphs  so  kj  assumes  only  discrete  values  spaced  uniformly  by  n/L  starting  at 
kz=0.  Using  the  same  dispersion  relation  for  the  cavity  as  for  the  infinite  guide,  the 
resonant  frequency  for  a  particular  mode  can  be  found  through  this  functional  relation  of  k^ 
and  to.  Therefore,  the  formulation  of  the  infinite  waveguide  can  be  experimentally  tested 
and  verified  in  the  laboratory  with  a  resonant  cavity!  The  major  benefit  of  a  resonant  cavity 
is  that  a  bounded  and  controlled  environment  is  provided  for  more  accurate  laboratory 
measurements. 

5.2  Resonant  Cavity  Measurements 

The  resonant  cavity  used  in  the  laboratory  measurements  was  fabricated  from 
aluminum  with  the  interior  dimensions  given  in  Figure  5.1.  With  the  excitation  and 
measuring  probes  located  in  the  centers  of  the  circular  end  caps,  the  TMqx  modes  are 
excited  f32].  Although  the  TM^  mode  is  not  fundamental  to  the  guide,  these  azimuthally 
symmetric  modes  are  easy  to  excite  and  prove  very  stable  in  the  sense  that  small  radial 
variations  in  the  probe  dimensions  do  not  exci.e  other  spurious  modes.  The  inside  of  the 
cavity  is  approximately  3.7  in  in  diameter  and  5.2  in  in  length.  The  measurements 
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performed  in  the  laboratory  were  for  the  three-dielectric  case:  a  hose  filled  with  some 
material.  The  hose  was  suspended  concentrically  along  the  cavity  axis  with  two  styrofoam 
rings  of  approximately  1/4  to  1/2  in  thickness.  It  was  determined  by  measuring  a  number 
of  resonant  frequencies  for  the  cavity,  both  empty  and  with  the  styrofoam,  that  the  rings 
provided  less  than  a  0.1%  change  in  the  center  frequency  of  the  resonance,  fc,  implying  the 
styrofoam’s  dielectric  constant  is  near  that  of  air.  The  loss  due  to  the  rings  was  also  small, 
since  Af3dB  increased  only  slightly.  The  measurements  of  the  cavity  Q  (fj/Afi^)  were  a 
source  of  great  difficulty.  The  primary  reason  for  this  was  the  lack  of  repeatability  in  the 
Af3dB  measurements  from  one  setup  of  the  cavity  to  the  next  and  from  one  day  to  the  next. 
Since  with  this  cavity  method  the  fc  measurements  were  very  consistent,  the  dielectric 
losses,  which  can  be  derived  from  the  Q,  were  measured  by  a  different,  more  repeatable 
method  to  be  discussed  in  Section  5.4. 

The  excitation  and  measuring  probes  in  the  cavity  were  straight,  32  gauge  wires 
protruding  into  the  cavity  0.13  in  from  the  end  cap.  From  trial-and-enor  measurements, 
both  the  diameter  and  length  of  wire  proved  to  be  critical  in  two  respects.  One,  too  large  a 
diameter  or  length  will  perturb  the  geometry  of  the  system  excessively  and  large  deviations 
from  the  simple  end-plate  boundary  conditions  <  Ep=E*=0)  become  evident  in  the  smoothed 

resonant  peaks.  Two,  for  a  probe  too  short,  the  coupling  to  the  fields  becomes  too  small 
and  the  resonances  are  lost  in  the  noise. 

The  source  of  excitation  for  the  cavity  was  a  Wiltron  6637A  signal  sweep 
generator.  To  measure  the  resonances,  a  Wiltron  560A  scalar  network  analyzer  was 
connected  to  the  cavity  output  through  type-N  connectors.  Figure  5.2  shows  the  output 
from  the  network  analyzer  for  the  arrangement  of  an  empty  0.46  in  inner  diameter  cpvc 
hose  of  thickness  0.0825  in  with  the  frequency  swept  from  2  to  10  GHz.  The  resonant 
peaks  are  quite  prominent  above  the  noise  level  near  -65  dB.  These  resonances  are  all 
associated  with  azimuthally  symmetric  TM  modes.  The  first  peak,  the  TM^o  resonance 
near  2.3  GHz,  corresponds  to  the  cutoff  frequency  of  the  TM01  mode  in  an  infinite 
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waveguide  of  similar  cross  section.  The  loaded  Q’s  for  these  lower-order  resonances  range 
from  a  low  of  approximately  300  to  a  high  near  1 100.  The  next  higher  mode,  the  TM^  in 
the  infinite  guide,  is  first  seen  in  the  stepped  region  near  5  GHz,  and  the  TM^  resonances 
begin  near  8  GHz  in  the  next  stepped  region.  Justification  for  these  designations  will 
become  evident  shortly  when  theoretically  predicted  resonances  are  compared  with  these 
results. 

When  „  hose  is  filled  with  a  silicone  sealant  (Dow  Coming  732),  the 
resonances  in  Figure  5.3  result.  All  the  peaks  experience  a  significant  downward  shift  m 
frequency,  but  the  higher-order  TM^  and  TM^  resonances  undergo  the  greatest  amount. 
A  comparison  of  the  resonant  peaks  for  the  above  two  cases  pius  those  of  the  empty  cavity 
is  shown  in  Figure  5.4.  Although  the  permittivities  of  the  cpvc  and  silicone  sealant  are 
relatively  small  (near  3.4  and  2.4  F/m,  respectively)  there  is  a  large  downward  shift  in  tne 
peaks  as  shown  in  the  figure.  Figures  5.5  and  5.6  show  the  resonant  peaks  of  the  same 
cavity  arrangement  but  having  a  3/4  in  cpvc  hose  of  inner  diameter  0.697  in  and  0.09  in 
thickness,  both  empty  and  filled  with  the  sealant,  respectively.  The  TMq,  resonances  of 
the  empty  3/4  in  hose  case  occur  lower  in  frequency  than  do  the  corresponding  ones  of  the 
1/2  in  hose  case.  Curiously,  for  both  the  TMq2  and  TM^  modes,  the  peaks  are  higher  in 
frequency  for  the  3/4  in  hose  than  for  the  1/2  in  one.  When  the  hose  is  filled  with  the 
sealant,  however,  all  the  resonances  of  the  3/4  in  hose,  including  the  higher-order  modes, 
occur  at  much  lower  frequencies  than  the  corresponding  resonances  of  the  1/2  in  hose.  All 
graphs  of  the  resonant  responses  shown  here  confirm  that  the  density  of  the  resonant  peaks 
increases  with  increasing  frequency,  making  it  quite  difficult  to  distinguish  the  individual 
peaks  of  the  higher-order  resonances. 
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5.3  Theoretically  Predicted  Resonances 

Comparison  of  the  measured  cavity  resonant  frequencies  with  the  theoretically 
calculated  ones  is  a  straightforward  process.  As  discussed  in  Section  5.1,  the  kz-co  plot 
for  the  infinite  waveguide  is  constructed,  then  the  appropriate  resonant  frequencies  of  the 
cavity  are  obtained  by  the  discretization  of  the  axial  wavenumber  into  its  allowed 
eigenvalues.  This  process  assumes  known  constitutive  parameters.  In  reality,  the  e.  s 
were  unknown  for  the  materials  used  in  the  laboratory.  Tnese  parameters  can  be  measured 
of  course,  but  due  to  a  lack  of  proper  equipment,  this  was  not  possible.  One  way  to  get 
around  this  stumbling  block  and  find  an  approximate  is  to  assume  a  value,  and  if  many 
calculated  resonant  frequencies  can  be  made  to  match  laboratory  measured  values,  then  one 
may  assume,  with  some  degree  of  certainty,  that  the  £,.  is  indeed  the  correct  parameter.  One 
difficulty  in  this  is  that  if  the  materials  are  dispersive,  as  most  are  in  the  microwave  region, 
the  Ej  would  need  to  be  continuously  varied  in  many  small  frequency  bands  to  accurately 
model  the  behavior.  Since  these  bands  are  unknown  to  begin  with,  it  is  not  possible  to 
accurately  model  dispersive  materials  using  this  method.  Another  difficulty  is  that  the 
materials  must  be  relatively  loss-free  since  this  peak-matching  method  is  unique,  or  nearly 
so,  only  in  the  loss-free  case.  Therefore,  it  will  be  assumed  here  that  the  materials  are 
lossless  and  nondispersive. 

For  the  cpvc  hose,  Fink  and  Christiansen  [33]  give  a  value  of  3.4  F/m  for  the 
dielectric  constant  at  1  GHz  and  tan  5=0.006.  Using  this  value  for  the  dielectric  constant 
across  the  entire  1  to  10  GHz  band,  the  predicted  resonant  frequencies  were  calculated  from 
the  (3—0)  plots  for  both  the  empty  1/2  and  3/4  in  lossless  cpvc  hoses.  Excellent  agreement 
was  realized  considering  for  32  peaks  which  clearly  matched  the  theoretical  and 
experimental  results,  there  was  an  average  error  of  less  than  1%  and  a  maximum  error  of 
2. 1  %.  For  the  filled  hose,  no  data  for  the  permittivity  were  available  on  the  silicone  sealant 
although  Dow  Coming  has  stated  that  £^.=2.8  F/m  with  tan  5=0.0015  at  100  kHz  [34], 
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Using  the  heuristic  peak  matching  method,  a  value  of  2.4  F/m  for  the  dielectric  constant 
was  arrived  at  for  the  sealant.  The  TM01,  TM^  and  TM^  kj-o)  plots  for  the  1/2  and  3/4 
in  cpvc  pipes  filled  with  a  material  of  ^=2.4  F/m  are  given  in  Figure  5.7.  For  29  peaks 
which  matched,  there  was  again  an  average  error  of  less  than  1%  and  a  maximum  error  of 
2.1%  with  the  measured  values.  A  selected  sample  of  values  for  these  two  sets  of 
comparisons  is  listed  in  Table  5.1.  More  than  likely,  the  permittivities  are  not  exactly  2.4 
F/m  for  the  sealant  and  3.4  F/m  for  the  cpvc  since  it  was  required  they  be  lossless  and 
nondispersive.  However,  it  is  probable  that  the  permittivities  are  very  near  those  given 
above  since  this  additional  boundary  condition  based  on  the  logical  argument  of  matching 
many  resonances  was  met  exceptionally  well. 

Finally,  the  resonant  cavity  technique  inherently  assumes  a  low-loss  material(s) 
within  the  cavity.  Too  high  a  loss  will  not  allow  the  cavity  to  resonate  since  by  the  time  the 
wave  traverses  the  length  of  the  cavity  twice,  the  signal  is  too  weak  to  contribute  to  the 
resonance  effect.  Measurements  on  methanol,  tap  water  and  glycerin  (all  very  lossy)  gave 
the  same  results —  the  measured  output  was  completely  devoid  of  resonant  peaks.  (A 
typical  example  of  this  is  shown  in  Figure  5.8  for  the  1/2  in  cpvc  hose  filled  with  methanol 
in  the  cavity.)  Other  techniques  are  needed  to  measure  the  loss  effects.  One  method  that 
gave  surprisingly  good  results  was  the  finite  waveguide  excited  by  horn  antennas  to  be 
discussed  next. 

5.4  Finite  Open-Ended  Waveguide  Measurements 

To  overcome  the  inherent  shortcomings  of  the  cavity-resonator  technique  s  inability 
to  measure  attenuation  for  very  lossy  dielectrics,  a  finite-waveguide  technique  was 
employed.  A  small  section  of  guide  was  attached  to  a  panel  on  a  shielded  enclosure,  as 
depicted  in  Figure  5.9,  and  illuminated  by  a  horn  antenna.  (A  more  detailed  description  of 
this  type  of  general  experimental  arrangement  is  given  in  [35].)  The  wave  attenuation 
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provided  by  the  guide  can  be  measured  with  regard  to  a  set  of  reference  measurements 
performed  with  the  two  antennas  in  free  space.  This  is  a  total  attenuation  from  all  sources 
including  diffraction,  reflection  and  wall  and  dielectric  losses.  The  antenna  gains  and 
system  losses  are  assumed  removed  by  the  reference  subtraction.  To  approximate  the 
losses  due  only  to  the  loading  dielectrics,  a  new  reference  level  can  be  chosen,  namely,  the 
total  guide  setup  with  the  empty  hose.  Then  upon  subtraction  of  this  reference  level  from 
the  guide  measurements  with  the  filled  hose,  the  diffraction  and  reflection  losses  will  be 
subtracted  in  addition  to  the  system  and  antenna  losses.  This  is  an  approximation  since  the 
fields,  in  addition  to  reflection  from  the  inner  filling,  will  be  redistributed  in  the  guide  from 
the  empty  to  filled  hose  cases.  Hence,  the  losses  from  other  sources  will  be  modified. 
However,  if  these  other  sources  of  loss  are  considered  negligible  and  the  filling  material 
inside  the  hose  the  only  source,  of  loss,  then  only  this  material  can  contribute  to  the 
attenuation  once  diffraction  and  reflection  arc  subtracted.  To  reduce  the  reflection  losses 
from  the  inner  material  (which  are  not  accounted  for  in  this  referencing  subtraction 
method),  conical-shaped  styrofoam  plugs  were  used  to  help  reflect  the  incident  energy 
down  the  guide.  All  these  seemingly  incredulous  assumptions  prove  quite  reasonable  for 
these  finite-guide  measurements  since  for  one,  the  guide  is  only  6  in  long  so  the  wall  losses 
are  small  and  two,  a  Teflon  hose  (Ejs2.l-j0.0003 15  F/m  at  3  GHz  [22])  is  used  to  contain 
the  material  and  is  quite  lossless.  Two  separate  examples  of  these  finite  guide 
measurements  will  be  discussed,  both  of  which  use  methanol  (which  is  very  lossy)  as  the 
filling  inside  a  Teflon  tube. 

The  first  example  is  for  a  3/4  in  Teflon  hose  of  inner  radius  0.3733  in  and  thickness 
0. 1908  in  filled  with  methanol  inside  a  1 .5  in  diameter  aluminum  waveguide  of  inner  radius 
0.7425  in  and  length  6  in.  The  results  of  this  comparison  arc  shown  in  Figure  5.10.  The 
measured  values  given  are  the  results  of  four  measurement  cycles  to  verify  the  repeatability 
of  the  experiment  and  identify  any  experimental  anomalies  such  as  humidity  effects.  The 
calculated  attenuation  shown  in  the  same  figure  is  for  the  HEn  mode  of  a  perfect  guide  of 
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similar  cross  section  and  length.  As  discussed  earlier,  this  eigensolution  method,  used  for 
the  calculated  values,  indicates  only  those  modes  that  can  exist  and  not  those  that  actually 
do  in  a  physical  arrangement.  However,  these  two  plots  are  remarkably  similar  which 
strongly  suggests  the  existence  of  this  mode  in  the  guide.  The  TM0i  mode  also  has  an 
attenuation  versus  frequency  graph  similar  to  this  one,  but  azimuthally  symmetric  modes 
are  not  likely  to  be  excited  by  the  incident  fields  from  the  pyramidal  horn  antenna  used 
here.  Most  likely,  there  are  many  modes  excited  in  addition  to  the  HEn,  but  here  it 
appears  that  the  major  behavior  of  the  attenuation  is  dictated  by  the  HEn  modal  properties 
alone. 

As  was  the  case  with  other  materials,  the  permittivity  values  for  the  methanol  are 
primarily  unknown;  however,  von  Hippel  [22]  lists  these  three  values 

er=30.9-j2.472  F/m  at  0.3  GHz 

e,=23.9-j  15.296  F/m  at  3.0  GHz 

er=8.9-j7.209  F/m  at  10.0  GHz  . 

Linear  interpolation  was  used  between  these  known  value"  and  is  shown  pictorially  in 
Figure  5. 1 1.  The  use  of  interpolation  here  is  a  bit  presumptive  since  the  real  and  imaginary 
parts  of  the  permittivity  cannot  be  chosen  independently.  These  quantities  are  related 
through  the  Kramers-Kronig  equations  which  specify  their  unique  behavior  for  proper 
satisfaction  of  the  causality  principle  [36],  [37],  [38]. 

Even  with  this  series  of  rough  assumptions,  the  agreement  between  the  theoretically 
calculated  attenuation  for  the  methanol  and  the  measured  values  is  remarkable.  Both  predict 
an  almost  linear  increase  in  attenuation,  on  the  log  scale,  as  the  frequency  is  increased.  For 
frequency  values  less  than  4.2  GHz,  the  air-filled  reference  guide  is  in  cutoff  and  no 
comparison  can  be  made 
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The  second  example  is  with  a  4  in  diameter  guide  of  radius  2.005  in  filled  with  a  1 
in  diameter  Teflon  tube  of  inner  radius  0.5008  in  and  thickness  0.1908  in.  The  same 
measurement  referencing  method  is  used  here  as  in  the  previous  case.  Excellent  agreement 
is  evident  again,  with  the  results  given  in  Figure  5.12  for  methanol  in  the  Teflon  tube. 
Surprisingly,  as  the  frequency  is  increased  here,  the  attenuation  peaks  near  6  GHz  and 
quickly  decreases  to  very  small  values  near  10  GHz.  This  behavior  could  not  be  modelled 
with  only  one  mode,  but  when  the  effects  of  the  EHn  mode  were  included  with  those  of 
the  HEu,  the  results  were  very  close  to  the  measured  values.  This  inclusion  of  the  EHn 
mode  indicates  a  strong  modal  conversion  between  the  HEn  and  EHn  modes  with  the 
EHn  dominating  at  higher  frequencies.  In  the  midband  region  where  neither  mode 
dominates,  the  correspondence  of  measured  to  predicted  values  is  marginal.  This  is  to  be 
expected  since  the  theoretically  calculated  values  assume  a  dominance  of  one  mode  over  all 
others;  clearly,  this  is  not  satisfied  in  the  midband  region.  Also,  above  8  GHz  the 
attenuation  becomes  negative  which  is  not  to  be  expected  for  this  passive  system.  This 
behavior  is  attributed  to  the  out-of-band  responses  of  the  pyramidal  horn  antennas  which 
have  single-mode  operation  only  up  to  8  GHz. 

A  number  of  observations  can  be  made  concerning  these  finite-guide 
measurements.  One,  the  EHn  mode  was  needed  in  the  4  in  guide  measurements  to 
improve  correspondence  with  measured  values  but  not  needed  for  the  1.5  in  guide. 
Examination  of  this  mode  in  the  small  guide  showed  that  the  minimum  attenuation  afforded 
was  near  200  dB  in  the  1  to  10  GHz  band.  Clearly,  this  mode  would  not  noticeably  affect 
the  attenuation  measurements.  When  it  does  dominate  as  in  Figure  5.12,  the  attenuation 
trend  changes  considerably.  This  introduces  the  next  observation:  the  attenuation  has  a 
negative  slope  and  approaches  a  small  value  in  the  4-in  guide  case  at  higher  frequencies. 
This  type  of  behavior  is  curious  from  the  standpoint  that  as  the  frequency  increases,  the 
attenuation  from  the  dielectric  losses  is  decreasing.  A  question  arises  as  to  whether  this  is 
reasonable  and  if  it  is  unique  to  this  physical  arrangement.  The  attenuation  as  a  function  of 
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frequency  for  the  homogeneous  waveguide  of  varying  radius  filled  with  methanol  was 
examined  using  the  linearly  interpolated  ^  values  for  methanol  as  shown  in  Figure  5.1 1.  It 
was  observed  that  attenuation  relationships  with  frequency  for  the  TEjj  and  TMn  modes 
assume  similar  shapes  as  in  the  inhomogeneous  guide  by  appropriate  choices  of  the  radius 
(not  necessarily  corresponding  to  either  guide),  but  at  much  higher  attenuation  (400-600 
dB).  However,  the  negative  slope  behavior  in  the  high  band  of  the  inhomogeneous  guide 
could  not  be  reproduced  with  the  homogeneous  guide —  the  attenuation  is  always 
increasing  in  this  band  or  asymptotically  so.  Although  the  conduction  losses  are  decreasing 
in  the  methanol  with  increasing  frequency,  in  the  homogeneous  guide  the  losses  remain  on 
the  rise,  while  in  the  inhomogeneous  4-in  guide  they  decrease.  This  decreasing  attenuation 
with  increasing  frequency  appears  to  be  unique  to  this  inhomogeneous  waveguide 


arrangement. 
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Measured  resonances  of  ihe  cavity  having  the  0.46  in  inner  diameter 
cpvc  pipe  of  thickness  0.0825  in  filled  with  the  silicone  sealant. 


Measured  cavity  resonances  when  filled  with  an  empty  cpvc 
pipe  of  inner  diameter  0.697  in  and  thickness  0.09  in. 
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Measured  resonances  of  the  cavity  having  the  0.697  in  inner  diameter 
cpvc  pipe  of  thickness  0.09  in  filled  with  the  silicone  sealant. 
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Figure  5.9.  Experimental  arrangement  for  the  finite-length  waveguide 
measurements  on  the  hose-filled  6  in  guide. 
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Figure  5. 10.  Comparison  of  the  theoretical  and  measured  attenuations  in  the  three- 
dielectric  guide  of  diameter  1 .485  in  containing  a  Teflon  lube  of  radius 
0.3733  in  and  thickness  0.1908  in  filled  with  methanol. 


Comparison  of  the  theoretical  and  measured  attenuations  in  the  three 
dielectric  guide  of  diameter  4.01  in  containing  a  Teflon  tube  of  radius 
0.5008  in  and  thickness  0.1908  in  filled  with  methanol. 


CHAPTER  6 

CONCLUSIONS  AND  SUGGESTIONS 
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FOR  FURTHER  STUDY 

In  this  thesis,  solutions  to  the  radially  inhomogeneous  infinite  circular  waveguide 
have  been  obtained  numerically  from  an  exact  theoretical  development  beginning  with  first 
principles.  These  inhomogeneous  waveguides  are  layered  with  circular,  concentric  annuli 
having,  perhaps,  complex  constitutive  parameters.  The  modes  that  were  found  to  exist  in 
these  waveguides  are  hybrid,  meaning  that  they  have  both  axial  E-  and  H-fields.  In 
general,  no  TE  and  TM  modes  are  possible  in  the  inhomogeneous  guide  except  for  modes 
with  n=0  and  all  modes  for  which  kz=0.  The  modal  derignation  scheme  and  modal 
nomenclature  for  these  hybrid  modes  have  not  yet  been  standardized  in  the  literature.  The 
modal  designation  technique  used  in  this  work  takes  the  homogeneous  guide  as  a  limiting 
case  such  that  as  the  constitutive  parameters  of  the  loading  dielectrics  approach  those  of  free 
space  or  the  dielectric  radii  are  increased  or  decreased  until  the  guide  becomes 
homogeneous,  the  hybrid  modes  approach  those  of  the  homogeneous  guide  as 

HEnm— >TEjUn  ,  EHnjn—^TMnm 

By  Waldron’s  Correspondence  Idea  [19],  these  hybrid  modes  have  a  1:1  correlation  with 
those  of  the  homogeneous  pec  guide.  Using  Waldron's  method  in  conjunction  with 
Snitzer's  scheme  of  modal  nomenclature  (the  sign  of  C11/C13),  a  unique  method  of 
defining  the  proper  modal  designation  is  obtained  which  is  both  useful  and  necessary  in  the 
tracing  process  for  kz.  This  tracing  process  can  be  complicated  by  such  occurrences  as 
mode  traces  intersecting  other  mode  traces  of  similar  azimuthal  variation  and  by  the 
development  of  backward-wave  regions  in  the  traces  when  the  permittivities  become 
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sufficiently  large  in  the  loading  dielectrics  relative  to  the  radius  of  the  dielectrics  and  the 
waveguide  wall. 

It  was  discovered  by  applying  this  method  of  modal  designation  and  nomenclature 
that  a  new  scheme  of  naming  the  backward-wave  modes  was  needed  since  older  methods 
led  to  nonphysical  results  such  as  unboundedness  in  the  waves  and  net  power  flow  directed 
toward  the  source  instead  of  away  in  the  backward-wave  regions,  as  examples.  This  new 
scheme  becomes  further  justified  after  a  small  amount  of  loss  is  introduced  into  the  loading 
dielectrics.  When  this  is  done,  the  complex  modes  are  found  to  be  coupled  forward  and 
backward  waves,  as  expected  from  the  conservation  of  energy  arguments.  After  increasing 
a/b  (or  some  other  physical  parameter)  to  a  sufficiently  large  value,  these  complex  modes 
decouple  into  forward  and  backward  waves.  As  to  whether  either  of  the  forward  and/or 
backward  waves  propagate,  additional  factors  must  be  considered,  such  as  the  relative  sizes 
of  the  Er,  values  and  relative  dielectric  radii.  This  idea  of  introducing  a  small  amount  of  loss 
to  separate  the  roots  is  a  very  useful  tool.  Here,  however,  it  complicates  matters 
considerably  since  the  complex  determinant  of  a  matrix  filled  with  Bessel  functions  of 
complex  arguments  must  be  computed. 

In  addition  to  using  dielectric  losses  as  tools  for  separating  roots,  these  complex  - 
constitutive  parameters  also  provide  attenuation  of  the  wave  which  is  especially  important 
in  shielding  applications.  In  this  study,  only  the  attenuation  provided  by  dielectric  losses 
was  examined  since  this  was  anticipated  to  be  the  major  source  of  attenuation  in  the 
inhomogeneous  guide  above  cutoff.  There  was  a  heavy  concentration  on  the  attenuation 
provided  by  the  lossy  three-dielectric  guide  since  this  is  a  good  model  for  the  hose  through 
the  waveguide-below-cutoff  penetration.  The  results  of  this  investigation  showed  that 
contrary  to  homogeneous  guides,  there  are  no  simple  rules-of-thumb  for  computing  the 
amount  of  attenuation  the  wave  suffers  in  relation  to  either  the  amount  of  dielectric  losses  or 
the  relative  sizes  of  the  dielectric  radii.  It  was  concluded  that  each  C3se  must  be  analyzed 
separately  by  tracing  out  the  axial  wavenumber  as  a  function  of  relative  dielectric  radius  or 
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frequency,  as  examples.  It  was  also  shown  that  the  redistribution  of  power  flow  in  the 
dielectrics  can  vary  the  attenuation  in  counterintuitive  ways.  For  example,  in  one  case, 
even  thougn  the  tosses  were  increased  considerably  in  the  inner  region  of  a  three-region 
waveguide,  the  attenuation  was  shown  to  decrease  appreciably. 

As  a  verification  of  the  theoretical  formulation  and  the  numerical  solutions,  two 
types  of  experimental  measurements  were  performed.  These  measurements  were  separated 
into  those  for  near  lossless  dielectrics  where  a  resonant  cavity  technique  was  employed  and 
those  for  lossy  materials  where  the  finite  guide  arrangement  for  the  waveguide-below- 
cutoff  penetration  was  used.  Both  gave  excellent  agreement  with  the  theoretically  predicted 
values.  In  the  finite  guide  arrangement,  a  number  of  speculative  assumptions  were  made 
concerning  the  reference  measurement  subtraction  from  the  complete,  loaded  waveguide 
arrangement  so  that  only  the  losses  from  the  loading  material  would  be  observed.  This  was 
necessary,  since  in  this  study,  only  the  characteristic  properties  of  the  waveguide  itself 
were  scrutinized  and  the  source  properties  were  not  considered. 

For  future  study  into  this  inhomogeneous  waveguide  topic,  an  important 
characteristic  of  this  penetration  which  deserves  considerable  attention  is  this  source 
coupling  to  the  waveguide.  By  considering  the  particular  source  of  EM  excitation  along 
with  the  waveguide  geometry,  the  total  penetration  arrangement  could  be  analyzed  such  that 
a  determination  of  the  actual  wave  attenuation,  including  diffraction  and  reflection  effects, 
could  be  calculated  in  a  deterministic  fashion  rather  than  identifying  excited  modes  and  the 
corresponding  attenuations  which  have  occurred.  This  would  allow  for  a  predictive  tool  in 
the  analysis  of  the  shielding  effectiveness  for  waveguide  penetrations. 

As  additional  study  into  this  phenomenon,  actual  measured  values  for  the 
permittivities  of  the  loading  dielectrics  could  be  obtained,  and  using  these  new  values  as 
input  to  the  tracing  programs,  more  accurate  comparisons  could  be  performed  with  the 
measured  field  quantities.  The  additional  boundary  condition  of  matching  many  resonant 
frequencies  applied  in  this  study  is  logically  consistent;  however,  using  actual  measured 
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permittivity  values  is  more  appealing  since  the  losses  present  in  the  loading  dielectrics  will 
have  some  effect  on  the  resonant  frequencies.  In  addition,  materials  with  dispersive 
constitutive  parameters  could  be  used  in  both  measurement  arrangements  with  a  higher 
degree  of  accuracy. 

Finally,  a  new  precision  resonant  cavity  fabricated  from  brass,  rather  that 
aluminum,  would  allow  for  more  accurate  resonant  frequency  measurements  in  addition  to 
repeatable,  accurate  Q  measurements.  From  the  Q,  another  measure  of  the  dielectric  losses 
could  be  correlated  with  the  theoretical  calculations. 


APPENDIX 


MDCW  COMPUTER  PROGRAM  LISTING 


The  program  MDCW  is  a  user  interactive  computer  program  for  finding  the  roots, 
the  resulting  field  coefficient  vector  and  the  transverse  fields  for  a  circular  waveguide  filled 
with  q  concentric,  layered  dielectrics.  Once  all  the  data  for  each  annulus  has  been  entered 
and  three  guesses  near  a  root  given,  the  program  fills  the  matrix  [A]  in  subroutine 
MTRXVAL.  Here,  by  virtue  of  Equations  (2.20)  through  (2.23),  values  are  added  row¬ 
wise  to  the  matrix.  These  equations  contain  the  Jn  and  Nn  Bessel  functions  which  are  not 
trivial  to  compute,  especially  for  complex  arguments.  Recursion  schemes  to  calculate  these 
Bessel  functions  suffer  from  excessive  error  corruption  if  upward  recursion  is  used  for 
more  than  a  few  orders.  Fortunately,  Donald  Amos  recently  published  a  package  of  higher 
transcendental  function  subroutines  which  give  incredibly  accurate  and  fast  results  without 
recursion  [39],  Instead,  u?  uses  a  number  of  series  approximations,  in  appropriate  ranges 
of  the  argument,  with  system  hardware  dependent  coefficients  to  maximize  the  available 
precision.  This  package  will  compute  Jn  and  N„  for  positive  n  and  complex  arguments  in 
addition  to  other  special  functions. 

The  computation  of  these  Bessel  functions  and  their  derivatives  is  carried  out  in 
subroutine  TOTBES.  Evaluation  of  the  determinant  of  the  filled  matrix  is  the  last  operation 
performed  in  subroutine  MTRXVAL.  Here,  as  discussed  earlier  in  Chapter  3,  LU 
decomposition  is  used  through  subroutines  DET1  and  LUDCMP.  It  is  the  value  of  the 
determinant  which  is  returned  to  the  main  program.  Once  the  determinant  becomes  small 
enough,  as  directed  by  the  root  finder,  signalling  the  location  of  a  root,  the  last  section  of 
the  main  body  of  MDCW  computes  the  transverse  fields  and  the  scalar  potentials  as  a 
function  of  the  waveguide  diameter  for  this  kz  root. 
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program  mdcwlO 


* 

* 

************************************************************************ 


Written  by  Keith  W,  Whites 

This  version  10  dated  March  16,  1988 


*  Program  --  Multiple  Dielectric  Circular  Waveguide  (mdcw) 

★ 

* 

* 

* 

* 

*  This  program  will  calculate  the  axial  wavenumber  for  a  circular  waveguide 

*  containing  q  concentric  dielectrics  having  complex  constitutive  parameters. 

*  Data  needed  by  the  program  for  these  calculations  are  frequency,  number  of 

*  dielectrics,  number  of  azimuthal  variations,  and  the  constitutive 

*  parameters  of  the  dielectrics  and  their  dimensions.  The  numerical  root 

*  finding  technique  used  is  Muller's  method  which  configures  three  points  for 

*  a  quadratic  fit  and  uses  a  root  of  that  quadratic  fit  to  move  towards 

*  the  zero  of  the  expression  to  be  evaluated  -  in  this  case  the  det[A], 


Output  Files:  fldmag.wvd  -  transverse  E  and  H  fields  (magnitudes,  phases, 

real  and  imaginary  parts)  as  a  function  of 
the  guide  radius. 

psie(m).wvd  -  scalar  electric  and  magnetic  potentials  as  a 
function  of  guide  radius. 

pwrin.wvd  -  input  file  for  "pwrflw3"  code  which  calculates 
the  total  transverse,  %  distribution  of  real, 
and  3-D  plot  of  power  flow. 

zkztra.wvd  -  trace  of  kz  (axial  wavenumber)  in  complex 
plane  as  the  root  finder  searches  for  the 
correct  root. 


*  Variables:  arg  -  argument  of  the  Bessel  functions 

*  besj  1  -  Bessel  function  of  the  first  kind  (J)  of  the 

*  inner  dielectric  of  the  two  under  evaluation 

*  besj2  -  J  of  the  outer  dielectric 

*  besn  1  -  Bessel  function  of  the  second  kind  (N)  of  the 

*  inner  dielectric  of  the  two  under  evaluation 

*  besn2  -  N  of  the  outer  dielectric 

*  bes**pl  -  Bessel  function  of  order  plus  one 

*  coeff  -  coefficients  [C]  for  the  field  expansions 

*  dbes**  -  derivative  of  the  Bessel  function  '**' 

*  dc  -  array  of  matrix  elements  for  which  kz  is  to  be 

*  determined  from 

*  freq  -  radian  frequency 

*  ipart  -  number  of  divisions  the  radius  is  divided  into  for 

*  the  field  and  scalar  wave  functions  plots 

*  ixl  -  matrix  row  pointer 

*  ix2  -  matrix  column  pointer 

*  nord  -  order  of  mode,  number  of  azimuthal  variations 

*  numdi  -  number  of  dielectrics 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

♦ 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

★ 

* 

★ 

* 

♦ 

* 

* 

* 

* 

★ 

* 

* 

* 

* 

★ 

* 

* 

* 

* 

* 
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psie(m)  -  scalar  electric  and  magnetic  potential  functions 
radii  -  radius  of  dielectrics  and  the  conductor  (m) 
sacc  -  accuracy  sought  for  k z  value 
umpl,ump2,ump3  -  determinant  values  for  three  guesses 
to  be  used  with  Muller’s  method 
uma,umb,umc  -  updating  information  for  Muller's  method 
zeps  -  permittivity 
zkro  -  radial  wavenumber 
zkz  -  axial  wavenumber 
zrou  -  permeability 


*  Subroutines: 


detl  -  calculates  the  determinant  of  a  double  precision 
complex  matrix  using  lu  decomposition 
lubksb  -  called  by  detl  to  perform  the  backsubstitution 
after  the  LU  decomposition  in  order  to  calculate 
the  coefficients  [C] 

mtrxval  -  calculates  the  elements  of  the  matrix  from 
the  boundary  conditions  for  the  multiple 
concentric  dielectric  circular  waveguide 
totbes  -  calculates  the  Bessel  functions  and  their 
derivatives  for  both  dielectrics  bordering 
an  interface 

zbesj(y)  -  called  by  totbes  to  calculate  the  Bessel 

functions  of  the  first  and  second  kinds  with 
complex  arguments 


************************************************************************ 


implicit  complex*  16(a-h,o-r,t-z) 
implicit  real*8(s) 

parameter  (pi=(3.14159265359d0,0.0d0),maxdie=5,maxcon=maxdie*4-2) 
parameter  (zeps0=(8.854d-12,0.0d0),zmu0=(12.56637062d-7,0.0d0)) 
parameter  (ipart=100,ibesnum=maxdie) 
complex*  16  zkro(maxdie)rradii(maxdie),zmu(maxdie), 

+  zeps(maxdie),dc(maxcon,maxcon),coeff(maxcon), 

+  efira(ipart),ef^hi(ipart),hfra(ipart),hfphi(ipart), 

+  psie(ipart),psieder(ipart),psim(ipart),psimder(ipart) 

complex*  16  worksp(maxcon),indx(maxcon) 
real*8  sjunkl(ibesnum),sjunk2(ibesnum) 
real*8  sbjr(ibesnum),$bji(ibesnum),sbyr(ibesnum),sbyi(ibesnurn) 
integer*4  ipivot(raaxcon) 
common/combe  s/nord, zkro /adii 
common/mtrx/freq,zmu, zeps, numdi, sacc.dc 


open(unit=  1 1  ,file=’a:zkztra.  wvd',status='unknown') 

open(unit=23,file='a:fldmag.wvd',status='unknown') 

open(unit=26,file=’a:pwrin.wvd',status='unknowri) 

open(unit=29,file=’a:psie.wvd',status=’unknown') 

open(unit=31,file='a:psim.wvd',status='unknown') 


:set  up 
:output  files 
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itcont=0 

mstp=0 

iflagd=l 

sguinc=O.OdO 

write(*,*)  'enter  frequency  (Hz)  and  the  number  of  dielectrics’ 

reaa(*,*)  sfreq,numdi 

freq=dcmplx(sfreq)*(2.0d0,0.0d0)*pi 

write(*,*)  'enter  order’ 

read(*,*)  nord 

write(*,*)  'enter  3  kz  guesses' 

write(*,*) '  (middle  guess, largest,smallest:recommended)' 

read(*,*)  zkzl 

szkzre=dreai(zkz  1 ) 

szkzim=dimag(zkz  1 ) 

write(ll,12)  szkzre,szkzim 

read(*,*)  zkz2 

szkzre=dreal(zkz2) 

szkzim=dkrag(zkz2) 

write(ll,12)  szkzre,szkzim 

read(*,*)  zkz3 

szkzre=dreal(zkz3) 

szkrim=dimag(zkz3) 

write(l  1,12)  szkzre,szkzim 

write(*,*)  'desired  accuracy?' 

read(*,*)  sacc 

* 

do  10  i=l,numdi 

write(*,*)  'enter  radius,  mu,  and  epsilon' 
read(*,*)  sradii,zmu(i),zeps(i) 
radii(i)=dcmplx(sradii) 
zmu(i)=zmu(i)*zmuO 
zeps(i)=zeps(i)*zepsO 
10  continue 

* 

write(*,*)  'Calculate  Field  Magnitudes?  (1/0)' 
read(*,*)  mageal 

* 

*  with  3  kz  guesses,  call  subrouune  mtrxval  to  calculate  matrix  elements 

*  and  the  value  of  the  resulting  determinant 

* 

call  mtrxval(zkzl,det,iflagd) 
umpl=det 

call  mtrxval(zkz2,det,iflagd) 
ump2=det 

call  mtrxval(zkz3,det,iflagd) 
ump3=det 

★ 

*  the  100  goto  loop  utilizes  Muller’s  method  to  determine  the  kz  value 

* 

100  itcont=itcont+l 

umq=(zkzl -zkz2)/(zkz2-zxz3) 

uma=umq*umpl-umq*((1.0d0,0.0d0)+umq)*ump2+umq*umq*ump3 
umb=((2.0d0,0.0d0)*umq+(  l  .0d0,0.0d0))*ump  1  -(( 1 .0d0,0.0d0)+umq) 
+  *((1.0d0,0.0d0)+umq)*ump2+umq*umq*ump3 


umc=((  1  .0d0,0.0d0)+umq)*ump  1 

* 

uradi=sqrt(umb*umb-(4.0d0,0.0d0)*uma*umc) 

* 

if  (abs(umb+uradi).gt.abs(umb-uradi))  then 
uradi=umb+uradi 
else 

uradi=umb-uradi 

endif 

* 

ump3=ump2 

ump2=umpl 

zkz3=zkz2 

zkz2=zkzl 

* 

*  calculate  a  new  point  and  update  last  two  points  for  quadratic  fit 

* 

zkz  1  =zkz2-  (zkz2- zkz3  )*  (2 .0d0,0 .  OdO)  *  umc/uradi 

szkzre=dreal(zkz  1 ) 

szkzjm=dimag(zlcz  1 ) 

write(ll,l2)  szkzre.szkzim 

call  mtrxval(zkz  1  .detjflagd) 

umpl=det 

* 

*  if  kz  varies  less  than  the  specified  value,  then  end 

* 

if  ((cdabs(zkzl-zkz2).lt.sacc))  then 

* 

write(V)  'finished' 
write(*,*)  'kz=',zkzl 

write(*,*)  'total  number  of  iterations  =  ’,itcont 
write(*,*)  'determinant  of  matrix  =  ',det 

★ 

do  200  i=l,numdi 
if(abs(zkro(i)).lt.  1  .OdO)  then 
write(*,*)  'potential  error :  zkro  small' 
write(*,*)  'zkro(',i,')=',zkro(i) 
endif 

write(*,*)  ’zkro(,,i,')=\zkro(i) 

200  continue 

* 

write(*,*)  'finished' 
mstp=l 
endif 

★ 

if(mstp)  1 00, 1 00,300 

* 

300  continue 
* 

if(magcal)999,999,305 

* 

305  continue 


maxnum=4*numdi-2 
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call  mtrxval(zkzl,det,0) 

* 

do  310  i=l,(maxnum-l) 
coeff(i)=-dc(i+l,l) 

310  continue 

* 

do  315  i=l,maxnum 
do  317  j=l,maxnum 
dc(ij)=dc(i+l  j+1) 

317  continue 
315  continue 
* 

*  calculate  field  expansion  coefficients 

* 

call  ludr  -np(dc,maxnum- 1  ,maxcon,indx,d) 
call  lubksb(dc,maxnum-  l,maxcon,indx,coeff) 
write(*,*)  'the  coefficients  are:' 

* 

do  320  i=maxnum.2,-l 
cocff(i)=coeff(i-l) 
write(*,*)  'coeff(',i,’)',coeff(i) 

320  continue 
* 

coeff(  1  )=(  1  .OdO.O.OdC) 
write(*,*)  'coeff(  l)',coeff(l) 
guidsp=radii(numdi)/dcmplx(ipax  t) 


.-shift  matrix  up 
-.and  to  the  left 
tone  element 


tarbitrarily  define 


*  calculate  transverse  E  and  H  fields  knowing  the  kz  and  coefficients 

★ 

do  400  j=l  apart 
sguinc=sguinc+dreal(guidsp) 

+ 

if(sguinc.lLdreal(radii(l)))  then 
arg=zkro(  1  )*dctnplx(sguinc) 

★ 

call  zbesj(drcal(arg),dirnag(arg),dble(nord),  1 ,2.sbjr,sbji, 

+  nz,ierr) 

iffierr.ne.O)  write(*,*)  Uessel  Routine  error  #’,ierr 
call  zbesy(dreal(arg),dimag(arg),dble(nord),l,2,sbyr,sbyi, 

+  nz,sjunkl,sjunk2,icrr) 

if(ierr.ne.O)  write(*,*)  Bessel  Routine  error  #',ierr 
bj=dcmplx(sbjr(  1  ))+(0.0d0, 1 .0d0)*dcmplx(sbji(  1 )) 
bjpl=dcmplx(sbjr(2))+(0.0d0,1.0d0)',,dcmplx(sbji(2)) 
dbj=dcmplx(nord)*bj/arg-bjp  1 

* 

efra(j)=-zkzl/(ffeq*zeps(  1  ))*coeff(  1  )*zkro(  1  )*dbj 
+  -dcmplx(nord)/dcmplx(sguinc)*coeff(2)*bj 

hfra(j)=-zkzl/(ffeq*zmu(l))*coeff(2)*zkro(l)*dbj 
+  -dcmplx(nord)/dcmplx(sguinc)*coeff(  1  )*bj 

efphi(j)=zkz  1  *dcmplx(nord)/(fieq*zeps(  1  )* 

+  dcmplx(sguir.c))*coeff(l)*bj+zkro(l)*coeff(2)*dbj 

hfphi(j)=-zkz  1  *dcmplx(nord)/(freq*zmu(  1  )* 

+  dcmplx(sguinc))*coeff(2)*bj-zkro(  1  )*coeff(  1  )*dbj 
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else 

* 

do  410  i=2,numdi 

* 

if(sguinc.ge.dreal(radii(i-l)).and.sguinc.le. 

+  dreal(radii(i)))  then 

img=i 
endif 

♦ 

410  continue 

* 

arg=zkro(  im  g)  *dcmplx(sguinc ) 

call  zbesj(dreal(arg),dimag(arg),dble(nord),l^,sbjr,sbji, 

+  nzden-) 

if(ierr.ne.O)  write(*,*)  'Bessel  Routine  error  #’,ierr 
call  zbesy(drcal(arg))dimag(arg),dble(nord),  1 ,2,sbyr,sbyi, 
+  nz,sjunkl,sjunk2,ierr) 

if(ierr.ne.O)  write(*,*)  'Bessel  Routine  error  #',ierr 
oj=dcmplx(sbjr(  1  ))+(0.0d0, 1  .OdO)*dcmplx(sbji(  1 )) 
by=dcmplx(sbyr(  1  ))+(0.0d0, 1  .OdO)*dcmplx(sbyi(  1 )) 

bjpl=dcmplx(sbjr(2))+(0.0d0, 1 ,0d0)*dcmplx(sbji(2)) 
byp  l=dcmplx(sbyT(2))+(0.0d0, 1 .0d0)*dcmplx(sbyi(2)) 
dbj=dcmplx(nord)*bj/arg-bjp  1 
dby=dcmplx(nord)*by/arg-byp  1 

ii=4*img-5 

e£ra(j)=-zkzl/(£req*zeps(img))*('coeff(ii)* 

+  zkro(img)*dbj+coeff(ii+l)*zktxj<img)*dby) 

cfra(j)=efra(j)-dcmpbc(nord)/dcmplx(sguinc)*(coeff(ii+2)* 
+  bj+coeff(ii+3)*by) 

hfra(j)=-zkzl/(freq*zmu(OTg))*(coeff(ii+2)* 

+  zkro(irng)*dbj+cocff(ii+3)*zkro(irng)*dby) 

hfira(j)=hfra(j>dcmplx(nord)/dcmplx(sguinc)*(coeff(ii)* 

+  bj+coeff(ii+l)*by) 

efphi(j)=zkz  1  *dcmplx(noid)/(freq*zeps(img)* 

+  dcmplx(sguinc))*(coeff(ii)*bj+coeff(ii+l)*by) 

efphi(j)=efphi(j)+zkro(img)*(coeff(ii+2)*dbj+ 

+  coeff(ii+3)*dby) 

hfphi(j)=-zkz  1  *dcmplx(nord)/(freq*  zmu(img)* 

+  ....  dcmplx(sguinc))*(coeff(u+2)*bj^oeff(ii+3)*by) 

hfphi(j)=hfphi(j)-zbx)(img)*(coeff(u)*dbj+ 

+  coeff(ii+l)*dby) 


psirn(j)=coeff(ii)*bj+coeff(ii+ 1  )*by 
psimder(j)=coeff(ii)*dbj+coeff(ii+l)*dby 
psie(j)=coeff(ii+2)*bj+coeff(ii+3)*by 
psieder(j)=coeff(ii+2)*dbj+coeff(ii+3)*dby 
endif 


400  continue 

* 

sfdmax=0.0d0 

* 

*  write  out  the  E  and  H  fields 

* 

do  420  i=l,ipart 

sfdmax=tnax(abs(efra(i)),abs(hfra(i)),abs(efphi(i)), 

+  abs(hfphi(i)),sfdmax) 

420  continue 

* 

sguinc=0.0d0 

write(23,*) '  radius  E  ro  (mag)  H  ro  (mag)  E  phi 
+  (mag)  H  phi  (mag)' 

* 

do  430  i=l,ipart 
sguinc=sguinc-Klreal(guidsp) 

write(23,24)  sguinc,abs(efra(i))/sfdmax,abs(hfra(i))/sfdmax, 
+  abs(efphi(i))/sfdmax,abs(hfphi(i))/sfdmax 

430  continue 

* 

sguinc=0.0d0 

write(23,*) r  radius  E  rc  (real)  H  ro  (real)  E  phi 
+  (real)  H  phi(real)' 

* 

do  431  i=l,ipart 
sguinc=sguinc+dreal(guidsp) 

write(23,24)  sguinc,dreal(efra(i))/sfdmax,dreal(hfra(i))/ 

+  sfdmax,dreal(efphi(i))/sfdmax,dreal(hfphi(i))/ 

+  sfdmax 

43 1  continue 

* 

sguinc=0.0d0 

write(23,*) '  radius  E  ro  (imag)  H  ro  (iraag)  E  phi 
+  (imag)  H  phi(imag)' 

* 

do  432  i=l,ipart 
sguinc=sguinc+dreal(guidsp) 

write(23,24)  sguinc,dimag(efTa(i))/sfdmax,dimag(hfra(i))/ 

+  sfdmax,dimag(efphi(i))/sfdmax,ciimag(hfphi(i))/ 

+  sfdmax 

432  continue 

* 

write(26,27)  dreal(zkzl),dimag(zkzl),sfreq,nord,numdi 

* 

do  440  i=l,numdi 

write(26,28)  dreal(zeps(i)/zeps0),dreal(zmu(i)/zmu0)t 
+  dreal(radii(i)) 

write(26,28)  dimag(zeps(i)/zeps0),dimag(zmu(i)/zmu0) 

440  continue 


do  450  i=l,maxnum 

write(26,28)  dreal(coeff(i)),dimag(coeff(i)) 
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450  continue 

* 

sguinc=0.0d0 

* 

*  write  out  the  electric  and  magnetic  scalar  potentials 

* 

do  460  i=l  ,ipart 

* 

sguinc=sguinc+dreal(guidsp) 

write(29,24)  sguinc,dreal(psie(i)),dimag(psie(i)), 

+  drcal(psieder(i)),dimag(psieder(i)) 

write(3 1 ,24)  sguinc,drcal(psim(i)),dimag(psim(i)), 

+  dreal(psimder(i)),dimag(psimder(i)) 

* 

460  continue 
* 

12  format(el2.5,5x,el2.5) 

24  format(el2.5,4x,el2.5,4x,el2.5,4x,el2.5,4x,el2.5) 

27  format(el2.5,5x,el2.5,5x,el2.5,5x,i4,5x,i4) 

28  format(el2.5,5x,el2.5,5x,el2.5) 

999  close(31) 

close(29) 

close(26) 

close(23) 

close(ll) 

end 

* 

* 

* 

* _ _ _ 

* 

*  Subroutine  mtrxval 

* 

*  mtrxval  calculates  the  matrix  elements  for  the  multiple  dielectric 

*  circular  waveguide  main  program. 

* 

*  input:  zkz  -  axial  wavenumber  guess 

*  iflag  -  =  1  for  determinant  calculation 

*  =  0  for  no  determinant 

* 

*  output:  det  -  determinant  of  matrix  for  kz  guess 

★ 

* _ _ 

★ 

* 

* 


subroutine  mtrxval(zkz,det,iflagd) 
implicit  complex*  16(a-h,o-r,t-z) 
implicit  real*8(s) 

parameter  (spi=3. 14159265359, maxdie=5,maxcon=maxdie*4-2) 
parameter  (zeps0=(8.854d-12,0.0d0),zmu0=(12.56637062d-7,0.0d0)) 
parameter  (ibesnum=4) 

complex*  16  zkro(maxdie)jadii(maxdie),zmu(maxdie), 

+  zeps(maxdie),dc(maxcon,maxcon) 
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complex*  16  worksp(maxcon) 

real*8  sjunkl(ibesnum),sjunk2(ibesnum) 

real*8  sbjr(ibesnum),sbji(ibesnum),sbyr(ibesnum),sbyi(ibesnum) 

integer*4  ipivot(maxcon) 

common/combes/nord,zkro,radii 

common/mtrx/freq,zmu,zeps,numdi,sacc,dc 

* 

do  1020  i=l,maxcon 
do  1030  j=l,maxcon 
dc(ij)=(0.0d0,0.0d0) 

1030  continue 

1020  continue 
* 

zkro(l)=sqn(freq*freq*zmu(l)*zeps(l)-zkz*zkz) 

Lx  1=1 
ix2=l 

* 

do  1000  i=2,numdi 

zkro(i)=sqn(&eq*£rcq*zmu(i)*zeps(i>zkz*zkz)  :principal  value 

arg==zkro(i)*radii(i-l) 

* 

*  for  the  given  zkz,  the  Bessel's  functions  and  their  derivatives  are  found 

★ 

call  totbes(i,arg,besj  1  ,besj2,bcsnl  ,besn2,dbesj  l,dbesj2, 

+  dbesnl,dbesn2) 

* 

*  calculate  matrix  elements,  inner  dielectric  is  a  special  case  and  is  source 

*  of  the  if  statements,  ixl-row  pointer,  ix2-column  pointer,  i=outer  dielectric 

*  of  the  two  under  observation. 

* 

dc(ix  1  ,ix2)=zkro(i- 1  )*zkro(i-  l)*zeps(i)*besj  1 

* 

iffixl.ne.l)  then 

dc(ix  1 4x2+1  )=zkro(i- 1  )*zkxo(i- 1  )*zeps(i)*besn  1 
ix2=ix2+2 
endif 

* 

dc(ix  1  4x2+2)=- zkro(i)*zkro(i)*zeps(i- 1  )*besj2 
dc(ixl,ix2+3)=-zkro(i)*zkro(i)*zeps(i-l)*besn2 
ixl=ixl  +  I 

* 

if(ixl.ne.2)  then 
ix2=ix2-2 

dc(ix  1  ,ix2+2)=zkro(i- 1  )*zkro(i- 1  )*zmu(i)*besj  1 
dc(ix  1 4x2+3)=zkro(i- 1  )*zkro(i-  n*zmu(i)*besn  1 
dc(ix  1  ,ix2+6)=-zkro(i)*zkro(i)*zmu(i- 1  )*besj2 
dc(ix  1  ,ix2+7)=-zJcro(i)*zkro(i)*zmu(i- 1  )*besn2 
else 

dc(ix  l,ix2+l  )=zkro(i-  l)*zkro(i- 1  )*zmu(i)*besj  1 
dc(ix  l,ix2+4)=-zkro(i)*zkro(i)*zmu(i-  l)*besj2 
dc(ix  1 4x2+5  )=-zkro(i)  *zkro(i)*zmu(i- 1  )*besn2 
endif 


ix  1  =ix  1+1 


dc(ixl,ix2)=zkz*dcmpix(nord)/(frcq*zeps(i- 1  )*radii(i- 1  ))* 
besjl 

if(ixl.eq.3)  then 
ix2=ix2-l 
else 

dc(ix  1  ,ix2+ 1  )=zkz*dcmplx(nord)/(freq*  zeps(i- 1 ) 
,"radii(i-l))*besnl 

endif 

dc(ix  1  ,ix2+2)=zkro(i- 1  )*dbesj  1 

if(ixl.eq.3)  then 
ix2=ix2-l 
else 

dc(ix  1  ,ix2+3)=zkro(i- 1  )*dbesn  1 
endif 

dc('x  1  ,ix2+4)=-zkz*dcmplx(nord)/(freq*zeps(i)*radii(i- 1  ))* 
besj2 

dc(ix  1  ,ix2+5)=-zkz*dcmplx(nord)/(freq*zeps(i)*radii(i- 1  ))* 
besn2 

dc(ix  1  ,ix2+6)=-zkro(i)*dbesj2 
dc(ix  1  ,ix2+7)=-zkro(i)*dbesn2 

if(ixl.eq.3)  then 
ix2=ix2+2 
endif 

ixl=ixl+l 

dc(ix  1  ,ix2)=zkro(i- 1  )*dbesj  1 

if(ixl.eq.4)  then 
ix2=ix2- 1 
else 

dc(ix  1  ,ix2+l  )=zkro(i- 1  )*dbesn  1 
endif 

dc(ix  1  ,ix2+2)=zJcz*dcmpIx(nord)/(ficeq*zmu(i- 1  )*radii(i- 1  ))* 
besjl 

if(ix  1  .eq.4)  then 
ix2=ix2-l 
else 

dc(ixl  ,ix2+3)=zkz<,dcmplx(nord)/(freq*zmu(i- 1 ) 
*radii(i-l))*besnl 

endif 

dc(ix  1  ,ix2+4)^-zkro(i)*dbesj2 
dc(ix  1  ,ix2+5)=-zkro(i)*dbesn2 

dc(ix  1  ,ix2+6)=-zkz*dempbc(nord)/(ffeq*zmu(i)*radii(i- 1  ))* 
besj2 

dc(ixl,ix2+7)=-zkz*dcmpU(nord)/(ffeq*zmu(i)*radii(i-l))* 

besn2 


if  (ixl.eq.4)  then 
ix2=ix2+2 
endif 

* 

ixl=ixl+l 

* 

if(ixl.eq.5)  then 
ix2=ix2+2 
else 

ix2=ix2+4 

endif 

* 

1000  continue 
* 

arg=zkro(numdi)*radii(numdi) 

* 

*  for  last  two  rows  calculate  elements,  zbesj(y)  is  Bessel's  routine  with  no 

*  derviatives.  last  two  rows  due  to  PEC  guiding  structure. 

* 

call  zbesj(dreal(arg),dimag(arg),dble(nord),  1 ,2,sbjr,sbji, 

+  nz.ierr) 

if(ierr.ne.O)  write(V)  'Bessel  Routine  error  #',ierr 
call  zbesy(dreal(arg),dimag(arg),dble(nord),l,2,sbyr,sbyi, 

+  nz,sjunkl,sjunk2,ierr) 

if(ierr.ne.O)  wri te(V)  'Bessel  Routine  error  #\ierr 
besj2=dcmplx(sbjr(  l))+(0.0d0, 1  .0d0)*dcmplx(sbji(  1 )) 
besn2=dcmplx(sbyr(  1  ))+(0.0d0, 1  .0d0)*dcmplx(sbyi(  1)) 
besj2pl=dcmplx(sbjr(2))+(0.0d0,1.0d0)*dcmplx(sbji(2)) 
besn2pl=dcmplx(sbyr(2))+(0.0d0,1.0d0)*dcmplx(sbyi(2)) 
dbesj2=dcmplx(nord)*besj2/arg-besj2p  1 
dbesn2=dcmplx(nord)*besn2/arg-besn2p  1 
maxnum=numdi*4-2 
dc(ix  1  ,maxnum-3)=besj2 
dc(ix  1  vmaxnum-2)=besn2 
ixl=ixl+l 

dc(ix  1  .maxnutn- 1  )=dbesj2 
dc(ix  1  ,maxnum)=dbesn2 
nc=numdi*4-2 

* 

*  with  matrix  calculated,  now  call  determinant  routine. 

* 

if  (iflagd)  2020,2020,2030 

* 

2030  continue 
* 

call  detl(dc,nc,det) 

* 

2020  continue 
return 
end 


Subroutine  totbes 


totbes  calculates  the  Bessel  functions  of  the  first  and  second  kinds 
and  the  corresponding  derivatives  for  complex  arguments. 


input:  i  -  loop  counter  (outer  dielectric  of  the  two  under  scrutiny) 
arg  -  argument  of  Bessel  functions 

output:  besjl,besj2  -  Bessel  functions  of  the  first  kind 

1 - inner  dielectric  of  the  boundary 

2- outer  dielectric  of  the  boundary 
besnl,besn2  -  Bessel  functions  of  the  second  kind 
dbesjl,dbesj2  -  derivative  of  the  Bessel  function  of  the 

first  kind 

dbesnl,dbesn2  -  derivative  of  the  Bessel  function  of  the 
second  kind 


subroutine  totbes(i,arg,besj  1  ,besj2,besn  1  ,besn2,dbesj  1 , 

+  dbesj2,dbesnl,dbesn2) 

implicit  complex*  16(a-h,o-r,t-z) 
implicit  rcal*8(s) 

parameter  (pi=(3.14159265359d0,0.0d0),raaxdie=5,maxcon=maxdie*4-2) 

parameter  (ibesnum=4) 

complex*  16  zkro(maxdie)jadii(maxdie) 

real*8  sjunkl(ibesnum),sjunk2(ibesnum) 

real*8  sbjr(ibcsnum),sbji(ibesnum),sbyr(ibesnum),sbyi(ibesnum) 

common/combe  s/nord,zkro/adii 

with  the  argument  and  the  order,  call  the  Bessel  function  subroutine 

call  zbesj(dreal(arg),dimag(arg),dble(nord),  1 ,2,sbjr,sbji, 

+  nz,ierr) 

if(ieiT.ne.O)  write(*,*)  'Bessel  Routine  error  #',ierr 
call  zbesy(dreal(arg),dimag(arg),dble(nord),  1 ,2,sbyr,sbyi, 

+  nz,sjunkl,sjunk2,ierr) 

if(ierr.ne.O)  write(*,*)  'Bessel  Routine  error  #',ierr 
besj2=dcmplx(sbjr(  1  ))+(0.0d0, 1  ,OdO)*dcmplx(sbji(  1 )) 
besn2=dcraplx(sbyr(  1  ))+(0.0d0, 1  .OdO)*dcmplx(sbyi(  1 )) 
besj2p  1  =dcmplx(sbjr(2))+(0.0d0, 1  .OdO)*dcmplx(sbji(2>) 
besn2p  1  =dcmplx(sbyr(2))+(0.0d0, 1 .0d0)*dcmplx(sbyi(2)) 
dbesj2=dcmplx(nord)*besj2/arg-besj2pl 
dbesn2=dcmplx(nord)*besn2/arg-besn2p  1 


calculate  Bessel  functions  and  derivatives  for  the  inner  dielectric 
of  the  boundary  under  consideration 

arg=zkro(i- 1  )*radii(i- 1 ) 

call  zbesj(dreal(arg),dimag(arg),dble(nord),l,2,sbjr,sbji, 

+  nz,ierr) 

if(ierr.ne.O)  write(*,*)  'Bessel  Routine  error  #',ierr 
call  zbesy(dreal(arg),dimag(arg),dble(nord),  1 ,2,sbyr,sbyi, 
+  nz,sjunkl,sjunk2,ierr) 

if(ierr.ne.O)  write(*,*)  ’Bessel  Routine  error  #',ienr 
besj  1  =dcmplx(sbjr(  1  ))+(0.0d0, 1  .OdO)*dcmplx(sbji(  1 )) 
besn  l=dcmplx(sbyr(  1  ))+(0.0d0, 1  .OdO)*dcmplx(sbyi(  1 )) 
besj  lp  1  =dcmplx(sbjr(2))+(0.0d0, 1  .OdO)*dcmplx(sbji(2)) 
besn  1  p  1  =dcmplx(sbyr(2))+(0.(WC,l  .0d0)*dcmplx(sbyi(2)) 
dbesj  l=dcmplx(nord)*besj  1/arg-besj  lpl 
dbesn  1  =dcmplx(nord)*besn  1/arg-besn  1  p  1 

return 

end 


Subroutine  detl 

detl  calculates  the  determinant  of  a  double  precision  complex 
matrix  using  lu  decomposition,  reference  "Numerical  Recipes." 


input:  a  -  square  matrix  under  evaluation 
n  -  size  of  square  matrix 

output:  d  -  determinant  (complex,  double  precision) 


subroutine  detl(a,n,d) 
implicit  complex*  16(a-h,o-z) 

parameter(maxdie=5,maxcon=maxdie*4-2,np=maxcon) 
complex*  16  a(np,np),indx(np) 
call  ludcmp(a,n,np,indx,d) 

do  1 1  i=l,n 
d=d*a(i,i) 


return 

end 


Subroutine  ludcmp 

ludcmp  performs  the  LU  decomposition  of  the  matrix  [dc]  to  be  used  in 
conjunction  with  the  root-finding  routine  in  searching  for  kz,  or 
in  conjunction  with  lubksb  to  calculate  coefficients  [C]. 
reference  "Numerical  Recipes." 

input:  a  -  square  matrix  under  evaluation 
n  -  size  of  square  matrix 
np  -  dimension  of  matrix 

output:  indx  -  row  permutations  for  partial  pivoting 
d  -  even  or  odd  number  of  row  interchanges 


subroutine  !udcmp(a,n,np,indx,d) 

implicit  complex*  16(a-h,o-z) 

parameter(nmax=  100) 

complex*  16  a(np,np)dndx(np),w(nraax) 

tiny=(  1 ,0d-40, 1 .0d-40) 

d=(1.0d0,0.0d0) 

do  12  i=l,n 
aamax=(0.0d0,0.0d0) 

do  11  j=l,n 

if  (cdabs(a(ij)).gt.cdabs(aamax))  aamax=abs(a(i,j)) 

1 1  continue 

if  (aamax.eq.O.)  pause  'singular  matrix' 
w(i)=(  1 .0d0,0.0d0)/aamax 

12  continue 

do  19  j=l,n 

do  14  i=l J-l 
sum=a(ij) 

do  13  k=lj-l 
sum=sum-a(ijc)*a(kj) 

13  continue 

a(ij)=sum 

14  continue 
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aamax=(O.OdO,O.OdO) 

do  16  i=j,n 
sum=a(i,j) 

do  15  k=l  j-1 
sum=sum-a(ijc)*a(kj) 

15  continue 

a(ij)=sura 

dum=w(i)*abs(sum) 
if(cdabs(dum).ge.cdabs(aamax))  then 
imax=i 
aamax=dum 
endif 

16  continue 

if  (j.ne.imax)  then 

do  17  k=l,n 
dum=a(imax,k) 
a(imax,k)=a(jjt) 
a(j,k)=dum 

17  continue 

d=-d 

w(imax)=w(j) 

endif 

indx(j)=imax 
if(a(j,j).eq.O.)  a(j,j)=tiny 
if(j.ne.n)  then 
dum=(  1 .0d0,0.0d0)/a(j  j) 

do  18  i=j+l,n 
a(i,j)=a(i,j)*dum 

18  continue 

endif 

19  continue 

return 

end 


Subroutine lubksb 

lubksb  calculates  the  backsubstitution  in  conjuction  with  ludemp 
to  calculate  the  coefficients  [CJ.  reference  "Numerical  Recipes." 


input:  a  -  square  matrix  under  evaluation 
n  -  size  of  square  matrix 
np  -  dimension  of  matrix 
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* 

* 

* 


* 

*  output:  indx  -  row  permutations  for  partial  pivoting 

*  b  -  right-hand  side  of  the  matrix  equation  Ax=b 

* 

* _ 

* 

* 

* 

subroutine  lubksb(a,n,np,indx,b) 
implicit  complex*  16(a-h,o-z) 
complex*  16  a(np,np),indx(np),b(np) 
ii=0 

* 

do  600  i=l,n 
ll=indx(i) 
sum=b(ll) 
b(ll)=b(i) 

* 

if  (ii.ne.O)  then 

* 

do  610  j=ii,i-l 
sum=sum-a(i,j)*b(j) 

610  continue 

* 

else  if  (sum.ne.0.)  then 
ii=i 
endif 

* 

b(i)=sura 

* 

600  continue 

* 

do  620  i=n,l,-l 
sum=b(i) 
if  (i.lt.n)  then 

* 

do  630  j=i+l,n 
sum=sum-a(i  j)*b(j ) 

630  continue 

* 

endif 

b(i)=sum/a(i,i) 

* 

620  continue 

* 


return 

end 
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